Supplementary information II
constants, kinetics and simulation code for simultaneously-fed, cross-templating systems



===
METHOD RK4							; Runge-Kutta 4th order - variable step methods are too fiddley
		{simult_sys_cross.mmd}
		{Random-time, *simultaneous* Gaussian spikes for *all* reactants, chemical & templated synthesis, then cross-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 poly (U) ms}

STARTTIME = 0							; housekeeping entries
STOPTIME = 600
DT = 1e-3							; spikes enter in 0.01 day
Tolerance = 1e-8						; low!
DTout = 1  							; partial output slightly faster
	

		{===rates & other system constants===} 
kd_pA = 5.3e-5							;  rate of decay of pN in *da^-1* units
kd_poly = 2.5e-3						; e.g., poly(U),  from decays in "Efficient heritable.."
kd_pN = 1e-2
kd_ImpA = 0.8             					; decay of activated pN or pR, da-1; instability due to activating group
kd_ImpN = 0.8
kd_AppA = 2.4e-3						; da-1, RppN, could have unstable reactive group
kd_AppN = 1e-2
kd_NppN = 2e-2
kt_AppA = 650							; M^-2 da^-1 3rd order reaction at templated nt interface; from "Efficient..." ms
kt_AppN = 650							; by analogy to AppA
kt_NppN = 650							; by analogy to AppA	
kc_AppA = 2.1							; /M /da, 2nd order reaction at solution stack interface; from poly(U) and "Efficient..." ms			
kc_AppN = 2.1							; by analogy to AppA
kc_NppN = 2.1							; by analogy to AppA
	
alpha = 0.5							; init fraction ImpA/pA+ImpA
gamma = 0.5       						; init fraction ImpN/pN+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

		{=== poly spikes ===}
spikesize_poly = 0.2						; mean total poly spikes, molar
spikeSD_poly = spikesize_poly/2					; SD of spike size; size/2

		{=== pN spikes ===}
TOTspikesize_pN = 0.2						; mean total pN spikes, molar
spikesize_pN = (1-gamma)*TOTspikesize_pN
spikeSD_pN = spikesize_pN/2					; SD of spike size, set for std pool behavior = 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-31

init pA = 1e-29
init poly = 1e-28
init pN = 1e-27
init ImpA = 1e-36
init ImpN = 1e-38
init AppA = 1e-35
init AppN = 1e-32
init NppN = 1e-30
	
init spiketime = 0
init spikeslope_pA = 0
init spikeslope_poly = 0
init spikeslope_pN = 0

init gaussspikedist_pA = 0
init gaussspikedist_poly = 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 -> 11, 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
	
		{=== poly-related spikes ===}
limit gaussspikedist_poly >= 0							; spike size >= 0
next gaussspikedist_poly = if spiketime >= 1 then gaussspikedist_poly 		; conserve spike size during spike
		else
		normal (spikesize_poly, spikeSD_poly)				; the Madonna way to normal variation; note BM error in "normal (mean, SD)": SD, not Var!
						
next spikeslope_poly = if spiketime >= 1 then gaussspikedist_poly/(10*DT) 	; conserve incoming slope during spike 
		else 
		0								; gaussian spike, converted to slope, and held const.during 11 spiketime intervals

		{=== 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_AppN*ImpN) -pA*poly*(kt_AppA*ImpA+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_AppN*pN) -ImpA*poly*(kt_AppA*pA+kt_AppN*pN)  															; reactant addn in 10 equal steps, with concurrent decay to pA
		{===integrate poly +/- spike===}
d/dt (poly) =  spikeslope_poly  -kd_poly*poly  					; template in 10 equal steps, with decay

		{===integrate pN +/- spike===}
d/dt (pN) =  spikeslope_pN +kd_ImpN*ImpN -kd_pN*pN -pN*(kc_AppN*ImpA+kc_NppN*ImpN) -pN*poly*(kt_AppN*ImpA+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_NppN*pN) -ImpN*poly*(kt_AppN*pA+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*poly*pA*ImpA

d/dt (AppN) = -kd_AppN*AppN +kc_AppN*(pA*ImpN+ImpA*pN) +kt_AppN*poly*(pA*ImpN+ImpA*pN)

d/dt (NppN) =  -kd_NppN*NppN +kc_NppN*pN*ImpN +kt_NppN*poly*pN*ImpN

d/dt (chem) =  +kc_AppA*pA*ImpA +kc_AppN*(pA*ImpN+ImpA*pN) +kc_NppN*pN*ImpN 																							; 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 			; chem for all products containing reactive N

d/dt (temp) =  +kt_AppA*poly*pA*ImpA +kt_AppN*poly*(pA*ImpN+ImpA*pN) +kt_NppN*poly*pN*ImpN				
										; total synthesis from solution templated stacks using G/A bp	
d/dt (temp_AppN) = +kt_AppN*poly*(pA*ImpN+ImpA*pN) 				; temp for NAD-like AppN alone

d/dt (temp_allN) = +kt_AppN*poly*(pA*ImpN+ImpA*pN) +kt_NppN*poly*pN*ImpN 	; temp for all products containing reactive N
	
tot = chem + temp								; total synthesis of NppN by all routes

ftemp = temp/tot