Supplementary information IV constants, kinetics and simulation code for simultaneously-fed, para-templating systems === METHOD RK4 ; Runge-Kutta 4th order - variable step methods are too fiddley {simult_sys_para,mmd} {simult_para: Random-time, */simultaneous/* Gaussian spikes for all reactants, chemical & template synthesis, then para-templated ribodimer synthesis, and decays; MY, 9/2017} {incorporates low-concentration form of rate equation, with 2nd and 3rd order rates for free and template-bound synthesis; uses expt'l rates from ggaa ms} STARTTIME = 0 ; housekeeping entries STOPTIME = 100 DT = 1e-3 ; spikes enter in 0.01 day Tolerance = 1e-8 ; low! DTout = 0.1 ; partial output slightly faster {===rates & other system constants===} kd_pA = 5.3e-5 ; rate of decay of pN in *da^-1* units kd_pG = 5.3e-5 kd_pN = 1e-2 kd_ImpA = 0.8 ; decay of activated pN or pR, da-1; instability due to activating group kd_ImpG = 1.3 kd_ImpN = 0.8 kd_AppA = 2.4e-3 ; da-1, RppN, could have unstable reactive group kd_GppG = 2.4e-3 kd_AppG = 2.4e-3 kd_AppN = 1e-2 ; instability due to intrinsic reactivity of pN kd_GppN = 1e-2 kd_NppN = 2e-2 kt_AppA = 27 ; M^-2 da^-1 3rd order reaction at templated nt interface; from GppG ms kt_GppG = 0 ; zero rates mean no catalysis of G nt on A templates kt_AppG = 0 kt_AppN = 27 ; by analogy to AppA kt_GppN = 0 kt_NppN = 27 ; by analogy to AppA kc_AppA = 2.1 ; /M /da, 2nd order reaction at solution stack interface; from poly(U) and GppG ms kc_GppG = 2.4 kc_AppG = 2.25 ; by interpolation kc_AppN = 2.1 kc_GppN = 2.25 ; by interpolation kc_NppN = 2.1 alpha = 0.5 ; init fraction ImpA beta = 0.5 ; init fraction ImpG gamma = 0.5 ; init fraction ImpN SpikeInt = 10 ; mean interval between spikes in life units totspikes = stoptime/SpikeInt ; # spikes of react in all, start -> stoptime threshold = totspikes*stepsize/stoptime ; spike threshold for totspikes using 0 -> 1 interval; probably 1st cycle -> 0 {=== pA spikes ===} TOTspikesize_pA = 0.2 ; mean total pA spikes, molar spikesize_pA = (1-alpha)*TOTspikesize_pA spikeSD_pA = spikesize_pA/2 ; SD of spike size, set for std pool behavior = size/2 {=== pG spikes ===} TOTspikesize_pG = 0.2 ; mean total pG spikes, molar spikesize_pG = (1-beta)*TOTspikesize_pG spikeSD_pG = spikesize_pG/2 ; SD of spike size; size/2 {=== pN spikes ===} TOTspikesize_pN = 0.2 ; mean total pN spikes, molar nt phosphate spikesize_pN = (1-gamma)*TOTspikesize_pN spikeSD_pN = spikesize_pN/2 ; SD of spike size; size/2 {===init condx===} init chem = 1e-26 init chem_AppN = 1e-26 init chem_allN = 1e-26 init temp = 1e-31 init temp_AppN = 1e-31 init temp_allN = 1e-25 init pA = 1e-29 init pG = 1e-28 init pN = 1e-27 init ImpA = 1e-36 init ImpG = 1e-37 init ImpN = 1e-38 init AppA = 1e-35 init GppG = 1e-34 init AppG = 1e-33 init AppN = 1e-32 init GppN = 1e-31 init NppN = 1e-30 init spiketime = 0 init spikeslope_pA = 0 init spikeslope_pG = 0 init spikeslope_pN = 0 init gaussspikedist_pA = 0 init gaussspikedist_pG = 0 init gaussspikedist_pN = 0 {===spike generator ===} next spiketime = if random (0,1) < threshold then 1 ; reset spiketime to give a spike else if spiketime >= 1 AND spiketime <= 9 then (spiketime + 1) ; internal spike count advances else 0 ; spiketime counts 1 -> 10, returns to 0 {=== pA-related spikes ===} limit gaussspikedist_pA >= 0 ; spike size >= 0 next gaussspikedist_pA = if spiketime >= 1 then gaussspikedist_pA ; conserve spike size during spike else normal (spikesize_pA, spikeSD_pA) ; the Madonna way to normal variation; note BM error in "normal (mean, SD)": SD, not Var! next spikeslope_pA = if spiketime >= 1 then gaussspikedist_pA/(10*DT) ; conserve incoming slope during spike else 0 ; gaussian spike, converted to slope, and held const.during 11 spiketime intervals spikeslope_ImpA = (alpha/(1-alpha))*spikeslope_pA ; activated slope simply alpha x pN; activation => constant, alpha {=== pG-related spikes ===} limit gaussspikedist_pG >= 0 ; spike size >= 0 next gaussspikedist_pG = if spiketime >= 1 then gaussspikedist_pG ; conserve spike size during spike else normal (spikesize_pG, spikeSD_pG) ; the Madonna way to normal variation; note BM error in "normal (mean, SD)": SD, not Var! next spikeslope_pG = if spiketime >= 1 then gaussspikedist_pG/(10*DT) ; conserve incoming slope during spike else 0 ; gaussian spike, converted to slope, and held const.during 11 spiketime intervals spikeslope_ImpG = (beta/(1-beta))*spikeslope_pG ; activated fraction => constant, beta {=== pN-related spikes ===} limit gaussspikedist_pN >= 0 ; spike size >= 0 next gaussspikedist_pN = if spiketime >= 1 then gaussspikedist_pN ; conserve spike size during spike else normal (spikesize_pN, spikeSD_pN) ; the Madonna way to normal variation; note BM error in "normal (mean, SD)": SD, not Var! next spikeslope_pN = if spiketime >= 1 then gaussspikedist_pN/(10*DT) ; conserve incoming slope during spike else 0 ; gaussian spike, converted to slope, and held const.during 10 spiketime intervals spikeslope_ImpN = (gamma/(1-gamma))*spikeslope_pN ;activated fraction => constant, gamma {===integrate pA +/- spike===} d/dt (pA) = spikeslope_pA +kd_ImpA*ImpA -kd_pA*pA -pA*(kc_AppA*ImpA+kc_AppG*ImpG+kc_AppN*ImpN) -pA*GppG*(kt_AppA*ImpA+kt_AppG*ImpG+kt_AppN*ImpN) ; reactant in 10 equal steps, with decay, ImpA -> pA accomodated in the mix {===integrate ImpA +/- spike===} d/dt (ImpA) = spikeslope_ImpA -kd_ImpA*ImpA -ImpA*(kc_AppA*pA+kc_AppG*pG+kc_AppN*pN) -ImpA*GppG*(kt_AppA*pA+kt_AppG*pG+kt_AppN*pN) ; ; reactant addn in 10 equal steps, with concurrent decay to pA {===integrate pG +/- spike===} d/dt (pG) = spikeslope_pG +kd_ImpG*ImpG -kd_pG*pG -pG*(kc_AppG*ImpA+kc_GppG*ImpG+kc_GppN*ImpN) ; reactant in 10 equal steps, with decay, ImpG -> pG accomodated in the mix {===integrate ImpG +/- spike===} d/dt (ImpG) = spikeslope_ImpG -kd_ImpG*ImpG -ImpG*(kc_AppG*pA+kc_GppG*pG+kc_GppN*pN) ; ; reactant addn in 10 equal steps, with concurrent decay to pG {===integrate pN +/- spike===} d/dt (pN) = spikeslope_pN +kd_ImpN*ImpN -kd_pN*pN -pN*(kc_AppN*ImpA+kc_GppN*ImpG+kc_NppN*ImpN) -pN*GppG*(kt_AppN*ImpA+kt_GppN*ImpG+kt_NppN*ImpN) ; reactant in 10 equal steps, with decay, ImpN -> pN in the mix {===integrate ImpN +/- spike===} d/dt (ImpN) = spikeslope_ImpN -kd_ImpN*ImpN -ImpN*(kc_AppN*pN+kc_GppN*pG+kc_NppN*pN) -ImpN*GppG*(kt_AppN*pA+kt_GppN*pG+kt_NppN*pN) ; ; reactant addn in 10 equal steps, with concurrent decay tp pN d/dt (AppA) = -kd_AppA*AppA +kc_AppA*pA*ImpA +kt_AppA*GppG*pA*ImpA d/dt (AppN) = -kd_AppN*AppN +kc_AppN*(pA*ImpN+ImpA*pN) +kt_AppN*GppG*(pA*ImpN+ImpA*pN) d/dt (NppN) = -kd_NppN*NppN +kc_NppN*pN*ImpN +kt_NppN*GppG*pN*ImpN d/dt (GppN) = -kd_GppN*GppN +kc_GppN*(pG*ImpN+ImpG*pN) d/dt (GppG) = -kd_GppG*GppG +kc_GppG*pG*ImpG d/dt (AppG) = -kd_AppG*AppG +kc_AppG*(pA*ImpG+ImpA*pG) d/dt (chem) = +kc_AppA*pA*ImpA +kc_AppN*(pA*ImpN+ImpA*pN) +kc_NppN*pN*ImpN +kc_GppN*(pG*ImpN+ImpG*pN) +kc_GppG*pG*ImpG +kc_AppG*(pA*ImpG+ImpA*pG) ; total synthesis from solution chemical stacks d/dt (chem_AppN) = +kc_AppN*(pA*ImpN+ImpA*pN) ; chem for NAD-like AppN alone d/dt (chem_allN) = +kc_AppN*(pA*ImpN+ImpA*pN) +kc_NppN*pN*ImpN +kc_GppN*(pG*ImpN+ImpG*pN) ; chem for all products containing reactive N d/dt (temp) = +kt_AppA*GppG*pA*ImpA +kt_AppN*GppG*(pA*ImpN+ImpA*pN) +kt_NppN*GppG*pN*ImpN ; total synthesis from solution templated stacks using G/A bp d/dt (temp_AppN) = +kt_AppN*GppG*(pA*ImpN+ImpA*pN) ; temp for NAD-like AppN alone d/dt (temp_allN) = +kt_AppN*GppG*(pA*ImpN+ImpA*pN) +kt_NppN*GppG*pN*ImpN ; temp for all products containing reactive N SpA = pA + ImpA + 2*AppA + AppN + AppG ; conservation, all forms of A, add decay SpG = pG + ImpG + 2*GppG + AppG +GppN ; conservation, all forms of G, add decay SpN = pN + ImpN + 2*NppN + AppN + GppN ; conservation, all forms of N, except for decay tot = chem + temp ; total synthesis of NppN by all routes ftemp = temp/tot ; fraction total synthesis as templating