Message Boards

GROUPS:

Staff Picks Biological Sciences Chemistry Mathematica Add-Ons Graphics and Visualization Wolfram Language Modeling Packages Statistics and Probability Publication Materials

William Piñeros

Computational framework for random chiral chemical networks

William Piñeros, University of Luxembourg

Posted 20 days ago

630 Views

1 Reply

10 Total Likes

Follow this post

WOLFRAM MATERIALS for the 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
Full article in PDF

enter image description here

Reaction schematic in a solution of achiral elements $A_i$(gray) and chiral elements $R_j (S_j)$, where $R_j$(blue) is the chiral counterpart of $S_j$ (red).

For more info see DOI: https://doi.org/10.1038/s41467-022-29952-8

Computational framework for random chiral chemical networks
by William Piñeros

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

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.

Model Set Up

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;

In[]:=

(*Reactioncontrol*)(*-------------------------------------------------*)Mrxn=Nelem+10;(*M:#ofrxnsinsystem.ShouldbeMrxn>Nelemtoavoidunderdetermination*)binprob=0.3;(*

:Probabilityofmakingbimolreaction*)catalprob=0.5;(*

:Probabilityofreactionbeingcatalyzedbyanindependentelement*)Δ=0.02;(*Δ:maximumpossiblenoise%inthesystem*)

In[]:=

(*forcingparameters*)(*-------------------------------------------------*)eta=0.6;(*η:Approximatefractionofreactionstobedriven*)f=10;(*f:totalnumberofavailableforces*)sJ=0.2;(*s:maximumforcingstrength*)mu=0.25;(*μ:Fractionofnon-zeroterms*)

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;

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.

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):

(

);(

driv

);

(symbolic)20-21. Separate stoichiometric matrix for R and S respectively (symbolic)

(*!!!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

In[]:=

(*Choosemodeltosimulate.The'modelSymBreak'canbeloadedbelowasanexample*)dummyModel=randomModel;dummyModel=modelSymBreak;

In[]:=

(*basicrunsetup.Tmaxshouldbelongenoughtoensuresteadystatesolution.~2-3

x10

isfineformostmodels.*)(*---------------------------------------------------------------*)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]];

In[]:=

(*---------------------------------------------*)(*---------------Simulation--------------------*)(*---------------------------------------------*)

In[]:=

(*Defineitoprocesss*)proc=ItoProcess[{XdotArgsDriven,noiseMat},{Xvars,Xvars0},t];(*Nowdothestochasticsimulation~~~maytakeafewmins,dependingonyoursystemandtmax*)nsims=1;simsD=RandomFunction[proc,{tmin,tmax,dt},nsims];

Analysis

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.

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}},JoinedTrue,PlotStyle{Black,Thick},LabelStyle{Large},FrameTrue,FrameLabel{"t(s)","Concentration"},MaxPlotPointsmaxpts],{elem,1,Nnormal}];chiralfigs=Table[td=trajExtract[sims,ti,tf,tstep,elem,Nchiral];ListLogLinearPlot[td,PlotRange{{ymin,ymax}},JoinedTrue,PlotStyle{Automatic,Dashed,Thick},LabelStyle{Large},FrameTrue,FrameLabel{"t(s)","Concentration"},MaxPlotPointsmaxpts],{elem,Nnormal+1,Neff}];Show[normfigs,chiralfigs]

Out[]=

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).

(*Selecttimerangeandsamplingpointstoaverageover.Theseshouldonlyfallwithinthesteady-statevaluesforaccuracy.Notethatinputsareactuallyindexessotimeisscaledbydt.*)ti=-110^3/dt;tf=-1;dtm=30/dt;(*Nonchiral*)drnormal=Table[test2a=trajSingExtract[simsD,ti,tf,dtm,nelem];m1=Mean[test2a];s1=StandardDeviation[test2a];{m1,s1},{nelem,1,Nnormal}];(*Chiralaverages*)drchiral=Table[test2a=trajSingExtract[simsD,ti,tf,dtm,nelem];test2b=trajSingExtract[simsD,ti,tf,dtm,nelem+Nchiral];m1=Mean[test2a];m2=Mean[test2b];s1=StandardDeviation[test2a];s2=StandardDeviation[test2b];{m1,s1,m2,s2},{nelem,Nnormal+1,Neff}];(*--------eeforaverageconcentrations-------*)ee[R_,S_]:=(R-S)/(R+S)*100;eed=Table[ee[drchiral[[i,1]],drchiral[[i,3]]],{i,Nchiral}];(*--------------------------Displaygrid---------------------------------*)(*Non-chiralelements*)meanColumn=PrependTo[drnormal,{"mean

","std

"}];elemColumn=PrependTo[Range[Nnormal],"Achiral elem #"];vertDiv=Thread[Range[3]True];horzDiv=Thread[Range[Length@meanColumn+1]False];horzDiv[[1]][[2]]=True;horzDiv[[2]][[2]]=True;horzDiv[[-1]][[2]]=True;NumberForm[Grid[Transpose@{elemColumn,meanColumn[[;;,1]],meanColumn[[;;,2]]},FrameAll,Dividers{vertDiv,horzDiv}],3](*chiral*)meanColumn=PrependTo[drchiral,{"mean

","std

","mean

","std

"}];eeColumn=PrependTo[eed,"ee %"];elemColumn=PrependTo[Range[Nchiral],"Chiral elem #"];vertDiv=Thread[Range[7]True];horzDiv=Thread[Range[Length@eeColumn+1]False];horzDiv[[1]][[2]]=True;horzDiv[[2]][[2]]=True;horzDiv[[-1]][[2]]=True;NumberForm[Grid[Transpose@{elemColumn,meanColumn[[;;,1]],meanColumn[[;;,2]],meanColumn[[;;,3]],meanColumn[[;;,4]],eeColumn},FrameAll,Dividers{vertDiv,horzDiv}],3]

Out[]//NumberForm=

Achiral elem #	mean A i	std A i
1	0.857	0.000177
2	0.361	0.000123

Out[]//NumberForm=

Chiral elem #	mean R i	std R i	mean S i	std S i	ee %
1	1.05	0.000272	0.44	0.000177	41.1
2	1.3	0.000607	0.645	0.000266	33.8
3	1.83	0.000517	0.362	0.000223	66.9
4	1.31	0.000207	0.244	0.000189	68.5
5	0.689	0.000364	0.507	0.000345	15.2
6	5.77	0.0014	0.172	0.000179	94.2
7	0.795	0.000273	1.92	0.000895	-41.4
8	0.545	0.00108	0.629	0.0012	-7.15
9	0.199	0.000181	0.715	0.000189	-56.5
10	0.483	0.000243	0.537	0.000344	-5.33
11	2.77	0.000542	0.399	0.000322	74.8
12	1.15	0.000156	0.323	0.000157	56.1

Work rate plots

Corresponding work rate plots of the above simulation. This is defined for each system channel R and S (see manuscript) as:



∑

where

is the net reaction of formation and

is the respective force acting on the reaction.For racemic systems, work rate is same on both channels seen as overlapping blue and orange lines. Chiral symmetry breaking systems will have diverging work rates, seen as a bifurcation in the lines.

In[]:=

(*Sampleofpointsforwhichtheworkrateiscalculated*)ti=.1;tf=tmax;tstep=1;samppoints=Range[ti,tf,tstep];(*Getabsoluatevalueofworkrateexpressions*)dotQR=dummyModel[[7]];dotQS=dummyModel[[8]];absdotQR=absTerms[dotQR];absdotQS=absTerms[dotQS];(*Computeworkrategivenconcentrationvaluesfromsimulation*){qtR,qtS}=dissSampTime[simsD,Xvars,absdotQR,absdotQS,samppoints];tt1=1;tt2=-1;dtt=1;ListLogLinearPlot[{qtR[[tt1;;tt2;;dtt]],qtS[[tt1;;tt2;;dtt]]},MaxPlotPointsmaxpts,JoinedTrue,PlotRangeAll,PlotStyle{Thick},LabelStyle{Large},FrameTrue,FrameLabel{"t (s) ","OverscriptBox[$\(w$\), $.$] (kT/s)"}]

Out[]=

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.

In[]:=

(*---Loadsamplemodel.FormattingsameasrandomModelgeneration(seeabove)---*)modelName="chiral_symmetry_breaking_example_model.m"modelSymBreak=Import[modelName];(*---Establishsizeofmodel---*)Nnormal=modelSymBreak[[-3]];Nchiral=modelSymBreak[[-2]];Mrxn=modelSymBreak[[-1]];

Out[]=

chiral_symmetry_breaking_example_model.m

(*****************************************************************************)
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
(*****************************************************************************)

POSTED BY: William Piñeros

1 Reply

Sort By:

Moderation Team

Moderation Team, WOLFRAM

Posted 20 days ago

-- you have earned *Featured Contributor Badge* Your exceptional post has been selected for our editorial column *Staff Picks* http://wolfr.am/StaffPicks and Your Profile is now distinguished by a *Featured Contributor Badge* and is displayed on the Featured Contributor Board. Thank you!

POSTED BY: Moderation Team

Reply to this discussion

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Group Abstract

Feedback