CITE THIS NOTEBOOK: Computational framework for random chiral chemical networks by William Piñeros. Wolfram Community JUN 10 2022.
ORIGINAL ARTICLE: Piñeros, W.D., Tlusty, T., Spontaneous chiral symmetry breaking in a random driven chemical system. NATURE Communications 13, 2244 (2022). https://doi.org/10.1038/s41467-022-29952-8.
ORIGINAL ARTICLE: Piñeros, W.D., Tlusty, T., Spontaneous chiral symmetry breaking in a random driven chemical system. NATURE Communications 13, 2244 (2022). https://doi.org/10.1038/s41467-022-29952-8.
Article abstract: Living systems have evolved to efficiently consume available energy sources using an elaborate circuitry of chemical reactions which, puzzlingly, bear a strict restriction to asymmetric chiral configurations. While autocatalysis is known to promote such chiral symmetry breaking, whether a similar phenomenon may also be induced in a more general class of configurable chemical systems—via energy exploitation—is a sensible yet underappreciated possibility. This work examines this question within a model of randomly generated complex chemical networks. We show that chiral symmetry breaking may occur spontaneously and generically by harnessing energy sources from external environmental drives. Key to this transition are intrinsic fluctuations of achiral-to-chiral reactions and tight matching of system configurations to the environmental drives, which together amplify and sustain diverged enantiomer distributions. These asymmetric states emerge through steep energetic transitions from the corresponding symmetric states and sharply cluster as highly-dissipating states. The results thus demonstrate a generic mechanism in which energetic drives may give rise to homochirality in an otherwise totally symmetrical environment, and from an early-life perspective, might emerge as a competitive, energy-harvesting advantage.
Introduction
Introduction
This notebook demonstrates the chiral chemical framework introduced in W.D. Pineros and T. Tlusty, Nature Comm. 13, 2244 (2022). In particular, this notebook:
1) generates instances of a random chemical network under an equally random set of driving forces;
2) simulates the corresponding dynamics and
3) provides analysis of the resulting work rate, steady state concentrations and enantiomeric excess of each of the elements, if any.
1) generates instances of a random chemical network under an equally random set of driving forces;
2) simulates the corresponding dynamics and
3) provides analysis of the resulting work rate, steady state concentrations and enantiomeric excess of each of the elements, if any.
Model Set Up
Model Set Up
Model parameters
Model parameters
Parameters explained in comments where symbol in italics denotes paper notation.
In[]:=
(*Numberofchiralandachiralelements*)(*-------------------------------------------------*)Nnormal=2;(*:Achiralelements*)Nchiral=12;(*:aChiralelementsR/S(getsdoubledbelow)*)Nelem=Nnormal+2*Nchiral;(*N:Totalnumberofelements*)Neff=Nnormal+Nchiral;
N
a
N
c
In[]:=
(*Reactioncontrol*)(*-------------------------------------------------*)Mrxn=Nelem+10;(*M:#ofrxnsinsystem.ShouldbeMrxn>Nelemtoavoidunderdetermination*)binprob=0.3;(*:Probabilityofmakingbimolreaction*)catalprob=0.5;(*:Probabilityofreactionbeingcatalyzedbyanindependentelement*)Δ=0.02;(*Δ:maximumpossiblenoise%inthesystem*)
p
b
p
c
In[]:=
(*forcingparameters*)(*-------------------------------------------------*)eta=0.6;(*η:Approximatefractionofreactionstobedriven*)f=10;(*f:totalnumberofavailableforces*)sJ=0.2;(*s:maximumforcingstrength*)mu=0.25;(*μ:Fractionofnon-zeroterms*)
Fixed parameters (leave unchanged)
Fixed parameters (leave unchanged)
In[]:=
(*ConcentrationconstantVolume=1,andonaverage1mol/vol*)Concentration=Nelem;(*Reactionconstants*)(*-------------------------------------------------*)kbase={1.,1.10^-1,1.10^-2,1.10^-3};(**)ksize=Length[kbase];(*Globallistofvarstobeusedformodelmaking*)(*-------------------------------------------------*)varlist=Nnormal,Nchiral,Mrxn,eta,Δ, binprob,catalprob,kbase,f,mu,sJ;
k
α
Generate model
Generate model
In[]:=
(*Offloadedtoamathematicapackageforautomateduse.PACKAGEDIR=locationofmathematicapackageinlocaldirectory*)(*---------------------------------------------------------------*)packagedir=NotebookDirectory[];PrependTo[$Path,packagedir];(*Invokepackage*)Needs["ChiralModelDemo`"];
Model generation: every time this cell is evaluated a brand new instance of the model is generated.
In this manner, one could generate, say, a table(array) of different models for mass evaluations etc. Here we focus on generating just one random instance of the model to be simulated next.
Note that, on average, chiral symmetry breaking models are very rare! Thus, you should expect that most generated models will be racemic. With default parameters, the expected rate of chiral symmetry breaking models is about ~3/200 random models, so a little bit of luck might be required.
In this manner, one could generate, say, a table(array) of different models for mass evaluations etc. Here we focus on generating just one random instance of the model to be simulated next.
Note that, on average, chiral symmetry breaking models are very rare! Thus, you should expect that most generated models will be racemic. With default parameters, the expected rate of chiral symmetry breaking models is about ~3/200 random models, so a little bit of luck might be required.
Generate instance of random model given chosen parameters. Most the of list contents are provided for completeness, and can be ignored. The few quantities used below for analysis will be automatically invoked where relevant.The returned list contains the following information:1. Element variable name (symbolic)2. Noise matrix (array)3. Strictly non-driven model (i.e. without forced terms) (symbolic)4. Driven + non-driven model. This is the full model to be simulated. (symbolic)5-6. Total forcing of the system, ||, for R and S channel respectively (symbolic)7-8. Total work rate for both channels as elaborated in below cells. (symbolic)9. List of reactions which were randomly chosen for driving. (list)10-11. Net reaction flux/currents for both channels respectively (symbolic)12-19. Rate constant vectors (for R and S respectively): (symbolic)20-21. Separate stoichiometric matrix for R and S respectively (symbolic)
∑
i
F
i
(,);(,);
+
k
-
k
+
k
driv
-
k
driv
(*!!!Beware,everyinvocationof'getModelAll'overridespreviousgeneratedmodel'randomModel',unlesssavedtoalistofmodels,arraysetc*)(*---------------------------------------------------------------*)Print["You've generated a brand new model"]randomModel=getModelAll[varlist];
You've generated a brand new model
Model simulation
Model simulation
In[]:=
(*Choosemodeltosimulate.The'modelSymBreak'canbeloadedbelowasanexample*)dummyModel=randomModel;dummyModel=modelSymBreak;
In[]:=
(*basicrunsetup.Tmaxshouldbelongenoughtoensuresteadystatesolution.~2-3isfineformostmodels.*)(*---------------------------------------------------------------*)tmin=0;tmax=910^3;dt=0.005;(*Symbolicvariablenames*)Xvars=dummyModel[[1]];(*Assigninitialvalueofvariablesfromrandomuniformdistributionbutsuchthat=*)Xvars0=randConcentrationChiral[Nelem,Nnormal,Nchiral,Nelem];(*Noisecoefficientmatrixforitosimulationbelow*)noiseMat=dummyModel[[2]];(*differentialequationexpressions,withforcings*)XdotArgsDriven=dummyModel[[4]];
4
x10
R
0
S
0
In[]:=
(*---------------------------------------------*)(*---------------Simulation--------------------*)(*---------------------------------------------*)
In[]:=
(*Defineitoprocesss*)proc=ItoProcess[{XdotArgsDriven,noiseMat},{Xvars,Xvars0},t];(*Nowdothestochasticsimulation~~~maytakeafewmins,dependingonyoursystemandtmax*)nsims=1;simsD=RandomFunction[proc,{tmin,tmax,dt},nsims];
Analysis
Analysis
Concentration Plots
Concentration Plots
Plots all elements as a function of time. Below: black lines = achiral elements; blue lines = chiral R_i; orange dashed = chiral S_i.
Clearly, if blue and orange lines overlap then the overall system is racemic. Chiral symmetry breaking instances would be observed as non-overlapping blue and orange lines.
Clearly, if blue and orange lines overlap then the overall system is racemic. Chiral symmetry breaking instances would be observed as non-overlapping blue and orange lines.
In[]:=
(*--ConcentrationPlot--*)sims=simsD;(*inputsasindeces*)ti=2;(*>2,otherwisewillgetanoriginartifactfromeplottingmethod*)tf=tmax/dt;tstep=1;tstep/=dt;(*tstep=timestep(s)toplot;Warning,smallerstepstakelongertoplot*)maxpts=2000;ymin=0;ymax=6;normfigs=Table[td=trajSingExtract[sims,ti,tf,tstep,elem];ListLogLinearPlot[td,PlotRange{{ymin,ymax}},JoinedTrue,PlotStyle{Black,Thick},LabelStyle{Large},FrameTrue,FrameLabel{"t(s)","Concentration"},MaxPlotPointsmaxpts],{elem,1,Nnormal}];chiralfigs=Table[td=trajExtract[sims,ti,tf,tstep,elem,Nchiral];ListLogLinearPlot[td,PlotRange{{ymin,ymax}},JoinedTrue,PlotStyle{Automatic,Dashed,Thick},LabelStyle{Large},FrameTrue,FrameLabel{"t(s)","Concentration"},MaxPlotPointsmaxpts],{elem,Nnormal+1,Neff}];Show[normfigs,chiralfigs]
Individual element stats and enantiomeric excess
Individual element stats and enantiomeric excess
Generate tables showing mean, std for achiral, and enantiomeric excess. This latter includes a computation of the average enantiomeric excess defined as:
ee = (R-S)/(R+S)*100
where R and S are mean concentrations values. The range of possible ee is [-1,1] where 0 = racemic and >0 for non-racemic distributions, with ±1 indicating total enantiomer dominance of either R or S respectively.
In practice, for racemic systems ee ~0 with small deviations due to noise and ee>>0 for chiral symmetry breaking systems (above noise level).
ee = (R-S)/(R+S)*100
where R and S are mean concentrations values. The range of possible ee is [-1,1] where 0 = racemic and >0 for non-racemic distributions, with ±1 indicating total enantiomer dominance of either R or S respectively.
In practice, for racemic systems ee ~0 with small deviations due to noise and ee>>0 for chiral symmetry breaking systems (above noise level).
Work rate plots
Work rate plots
Sample chiral symmetry breaking model
Sample chiral symmetry breaking model
Because the probability of getting a symmetry breaking model from random chance is low, we include the representative sample model used in the manuscript. This can be loaded, simulated and compared to the same analysis as above.
(*****************************************************************************)
Now can go back to the previous cells and re-simulate and analyze using this model in the ‘Model Simulation’ cell. Simulation should go for at least tmax=1x10^4 to ensure symmetry breaking is observed
(*****************************************************************************)
Now can go back to the previous cells and re-simulate and analyze using this model in the ‘Model Simulation’ cell. Simulation should go for at least tmax=1x10^4 to ensure symmetry breaking is observed
(*****************************************************************************)