We formulate a 3D ODE for photosynthesis that explicitly tracks the concentrations of [R], [ATP], and [H+] (the proton concentration, which modulates jPQ. We use mass conservation laws to track the concentrations of [ADP], [Pi], and [PGA]:
Conservation of phosphate: C1=2∗[R]+[PGA]+[Pi]+[ATP] Exchange of ADP/ATP: C2=[ATP]+[ADP] Equation for photochemical quenching flux jPQ: jPQ=kPQ+kf+kNPQkPQI where I is the incident light flux and the terms in the numerator and denominator are the probabilities of photochemical quenching, flourescence, and NPQ. Thus kPQ=1−kNPQ−kf. kNPQ depends on [H+]: kNPQ=(1−kf)kNPQ∗[H+]+1kNPQ∗[H+] Thus k_{NPQ} will increase with [H+] and saturate at 1−kf.
ATP is synthesized via the ATP Synthase SU, which processes ADP, Pi, and protons to produce ATP: $$ j_{ATP-syn} = \frac{[E_{syn}]}{\frac{1}{k_{syn}} + \frac{3}{f([ADP])} + \frac{12}{f([H^+])} + \frac{3}{f([P_i])} - \frac{1}{f([ADP])/3+f([H^+])/12}
\frac{1}{f([ADP])/3 + f([P_i])/3} - \frac{1}{f([P_i])/3+f([H^+])/12}
\frac{1}{f([ADP])/3+f([P_i])/3+f([H^+])/12}} $$ Where $f(A)isthediffusionrateofsubstrateAwithbindingratek^+A:f(A) = k^+A[A].LATEXE{syn}istheconcentrationofATPSynthase,andk{syn}isitsrateofproduction.<br>RuBPisconsumedviaCarboxylationandOxygenationandisregeneratedviatheRuBPregenerationSU,whichalsoconsumesATPandregeneratesPiand[ADP](weassumeherethatalloftheorthophosphatethatisconsumedtoproduceATPisregeneratedalongwithRuBP.Wedonotyetexplicitlymodelstarchandsucrosesynthesis.)jcA=dtd[C]=1/kERC+1/(f(O)+f(C))+1/f(R)−1/(f(C)+f(O)+f(R))EtLATEXjoA=dtd[O]=1/kERO+1/(f(O)+f(C))+1/f(R)−1/(f(C)+f(O)+f(R))Etandj{PGA-prod} = 2j_{cA} + 1.5j_{oA}.∗(NotethatE_tandk_{ERO}likelyneedstoapproximatethereactionsoftheoxidativepathwayinadditiontoCarboxylation)∗<br>TheRuBPregenerationSUactsasasimplifiedSUofthemanyenzymaticreactionsthateventuallyconvertPGAbacktoRuBP.2molPGAand3molATPproduce1molRuBP.jRuBP−regen=1/jCC+3/f(ATP)+2/f(PGA)−(f(ATP)/3+f(PGA)/21ECCThuswehavethefollowingODEsystem:dtd[R]=jRuBP−regen−158jPGA−prod+jleakwherej_{leak}istheRuBPlosttooxygenation.dtd[ATP]=3jATP−syn−3jRuBP−regenLATEXdtd[H+]=3jPQ−12jATP−syn$
from scipy import integrate as spint from scipy.optimize import fsolve from sympy import Matrix import numpy as np import matplotlib.pyplot as plt import ipywidgets as widgets from ipywidgets import interact, interact_manual inv = np.linalg.inv norm = np.linalg.norm eye = np.identity
#Parameter values #from Farquhar 1980 #turnover numbers for carboxylation and oxygenation kERC = 2.5 kERO = .21*kERC #turnover rate for carboxylation (1/s) kc = 4.14*((10)**(-3)) #turnover rate for oxygenation (1/s) ko = 7.57*(10**(-6)) #RuBP binding rate kr = 2*(10**(-5)) #light intensity (umol) I = 1000 #total enzyme concentration (umol/carboxylation site) Et = 87.2 #CO2 concentration (ubar) C = 230 #O2 concentration (ubar) O = 210 * 10**(3) kprod = kc kE = kc kappac = 2/3 kappar = 1/3 ATPo = 200 Ro = 350 ADPo = 100 Po = 100 PGAo = 200 Ho = 100 #total ATP + ADP (umol) OP = ATPo + ADPo OPt = ATPo + Po #total RuBP + intermediates (umol) tot = 2*Ro + OPt + PGAo Ecalvin = 100 Esyn = 100 kcalvin = 2 ki = 2 kAC = 2 kATP = 2 jsyn = 2 kp = 2 knpq = .02 kf = .05 kh = 2
#ODE system def Hmodel(x, I = I, knpq = knpq, C=C, O=O, kr = kr, ATPo = ATPo, Ro = Ro, PGAo = PGAo, Po=Po, ADPo=ADPo): R, ATP, H = x[0], x[1], x[2] Hd = knpq*H fnpq = (1-kf)*(Hd/(1+Hd)) jpq = (1 - (kf + fnpq))*I #total ATP + ADP (umol) OP = ATPo + ADPo OPt = ATPo + Po #total RuBP + intermediates (umol) tot = 2*Ro + OPt + PGAo P = OPt - ATP INte = tot - (ATP + 2*R + P) ADP = OP - ATP fC = kc*C fO = ko*O fR = kr*R #new = """ jcA = Et/((1/kERC) + (1/(fC+fO)) + 1/fR - 1/(fC+fO+fR)) joA = Et/((1/kERO) + (1/(fC+fO)) + 1/fR - 1/(fC+fO+fR)) jPGAprod = 2*jcA + 1.5*joA #""" #ATP consumption jRuBPregen = Ecalvin/(1/kcalvin + 3/(kAC*ATP) + 2/(ki*INte) - 1/(ki*INte/2+kAC*ATP/3)) #Pi consumption/ATP synthesis fATP = kATP*ADP/3 fP = kp*P/3 fH = kh*H/12 jPcons = Esyn/(1/jsyn + 1/fATP + 1/fH + 1/fP - (1/(fH+fATP)+1/(fH+fP)+1/(fP+fATP)) + 1/(fH+fP+fATP)) #this is the amount of phosphate lost due to oxygenation. See below leak = 232533.333333333/(20000/(20000*R + 2.54190000000000) - 1/R - 15868.1301388725) + 174400.000000000/(20000/(20000*R + 2.54190000000000) - 1/R - 45963.3682341106) cons = [jPGAprod, jRuBPregen, jPcons, jpq, leak] S = [[-8/15, 1, 0, 0, -1/2], [0, -3, 3, 0, 0], [0, 0, -12, 3, 0]] DEs = np.matmul(S, cons) return DEs, jpq
vs = ["R", "ATP", "H"] dt = .001 steps = 10000
def findPGA(R, ATP, Ro=Ro, ADPo=ADPo, PGAo=PGAo, ATPo=ATPo, Po=Po): OPt = ATPo + Po P = OPt - ATP tot = 2*Ro + OPt + PGAo PGA = tot - (ATP + 2*R + P) return PGA
#numerical simulation @interact(dt = (.00001, .01, .00001), I=(0,2000,10), steps=(10, 50000, 1), knpq = (.001, 1, .001), kr = (.01, 20000, .01), Ho = (0,100, 1)) def simiter3(steps = steps, dt = dt, Ro = Ro, Ho=Ho, ATPo=ATPo, PGAo = PGAo, Po=Po, ADPo=ADPo, I=I, O=O, C=C, knpq=.5, kr = kr): x0 = [Ro, ATPo, Ho] a = len(x0) sol = np.zeros((steps, a)) sol[0] = x0 jpqlist = [] for i in range(steps-1): DE, jpq = Hmodel(sol[i], I=I, knpq=knpq, C=C, O=O, kr =kr, ATPo = ATPo, Ro = Ro, PGAo = PGAo, Po=Po, ADPo=ADPo) step = DE*dt sol[i+1] = np.add(sol[i], step) jpqlist.append(jpq) print(sol[-1]) check, jp = Hmodel(sol[-1], I=I, knpq=knpq, C=C, O=O, kr =kr, ATPo = ATPo, Ro = Ro, PGAo = PGAo, Po=Po, ADPo=ADPo) print(check, jp) PGAlist = [findPGA(x[0], x[1], Ro=Ro, ADPo=ADPo, PGAo=PGAo, ATPo=ATPo, Po=Po) for x in sol] for k in range(a): y = [x[k] for x in sol] plt.plot(range(steps), y, label="{}".format(vs[k])) R = sol[-1][0] fC = kc*C fO = ko*O fR = kr*R jcA = Et/((1/kERC) + (1/(fC+fO)) + 1/fR - 1/(fC+fO+fR)) print(jcA) plt.plot(range(steps), PGAlist, label="PGA") plt.plot(range(steps-1), jpqlist, label = "jpq") plt.legend() plt.show()
#Loss due to oxygenation per cycle def oxyloss(): R, ATP, H, P, INte, ADP = var("R, ATP, H, P, INte, ADP") Hd = knpq*H fnpq = (1-kf)*(Hd/(1+Hd)) jpq = (1 - (kf + fnpq))*I fC = kc*C fO = ko*O fR = kr*R #new = """ jcA = Et/((1/kERC) + (1/(fC+fO)) + 1/fR - 1/(fC+fO+fR)) joA = Et/((1/kERO) + (1/(fC+fO)) + 1/fR - 1/(fC+fO+fR)) jPGAprod = 2*jcA + 1.5*joA #""" #ATP consumption jRuBPregen = Ecalvin/(1/kcalvin + 3/(kAC*ATP) + 2/(ki*INte) - 1/(ki*INte/2+kAC*ATP/3)) #Pi consumption/ATP synthesis fATP = kATP*ADP/3 fP = kp*P/3 fH = kh*H/12 jPcons = Esyn/(1/jsyn + 1/fATP + 1/fH + 1/fP - (1/(fH+fATP)+1/(fH+fP)+1/(fP+fATP)) + 1/(fH+fP+fATP)) leak = 232533.333333333/(20000/(20000*R + 2.54190000000000) - 1/R - 15868.1301388725) + 174400.000000000/(20000/(20000*R + 2.54190000000000) - 1/R - 45963.3682341106) cons = [jPGAprod, jRuBPregen, jPcons, jpq, leak] S = [[-8/15, 1, 0, 0, -1/2], [0, -3, 3, 0, 0], [0, 0, -12, 3, 0], [0, 3, -3, 0, 0], [1, -2, 0, 0, 0], [0, 3, -3, 0, 0] ] DEs = np.matmul(S, cons) dC = 2*DEs[0] + DEs[1] + DEs[3] + DEs[4] #leakage is equal to difference between increase in 2*R and that of PGA. So should maybe see under what cases one is greater than another? #Well it's 2R + PGA so #RuBP has two phosphates. ATP gives one phosphate and obviously Pi has one. Pretty sure PGA has three though?? PGA has one. Two moles of PGA are made from one RuBP Rlist = np.linspace(0.1, 10, 100) dClist = [dC.subs(R==r) for r in Rlist] #leaklist = [leak.subs(R==r) for r in Rlist] #diff = dC - leak #difflist = [diff.subs(R==r) for r in Rlist] plt.plot(Rlist, dClist, label='dClist') #plt.plot(Rlist, leaklist, label='leaklist') #plt.plot(Rlist, difflist, label = 'difflist') plt.title("Leakage with correction factor") plt.xlabel("RuBP concentration") plt.ylabel('Leakage') plt.legend() plt.show() oxyloss()
#now that we've added the correction factor to add a small amount to RuBP, the leakage is on a scale of 10^-14, likely due to rounding errors def leak(R): leak = 232533.333333333/(20000/(20000*R + 2.54190000000000) - 1/R - 15868.1301388725) + 174400.000000000/(20000/(20000*R + 2.54190000000000) - 1/R - 45963.3682341106) return leak
#Jacobian for backward Euler def dF(x): R, ATP, H = var("R, ATP, H") Hd = knpq*H fnpq = (1-kf)*(Hd/(1+Hd)) jpq = (1 - (kf + fnpq))*I #total ATP + ADP (umol) OP = ATPo + ADPo OPt = ATPo + Po #total RuBP + intermediates (umol) tot = 2*Ro + OPt + PGAo P = OPt - ATP INte = tot - (ATP + 2*R + P) ADP = OP - ATP fC = kc*C fO = 2*ko*O/3 fR = kr*R #new = """ jcA = Et/((1/kERC) + (1/(fC+fO)) + 1/fR - 1/(fC+fO+fR)) joA = Et/((1/kERO) + (1/(fC+fO)) + 1/fR - 1/(fC+fO+fR)) jPGAprod = 2*jcA + 1.5*joA #""" #ATP consumption jRuBPregen = Ecalvin/(1/kcalvin + 3/(kAC*ATP) + 2/(ki*INte) - 1/(ki*INte/2+kAC*ATP/3)) #Pi consumption/ATP synthesis fATP = kATP*ADP/3 fP = kp*P/3 fH = kh*H/12 jPcons = Esyn/(1/jsyn + 1/fATP + 1/fH + 1/fP - (1/(fH+fATP)+1/(fH+fP)+1/(fP+fATP)) + 1/(fH+fP+fATP)) leak = 232533.333333333/(20000/(20000*R + 2.54190000000000) - 1/R - 15868.1301388725) + 174400.000000000/(20000/(20000*R + 2.54190000000000) - 1/R - 45963.3682341106) cons = [jPGAprod, jRuBPregen, jPcons, jpq, leak] S = [[-8/15, 1, 0, 0, -1/2], [0, -3, 3, 0, 0], [0, 0, -12, 16, 0]] DEs = np.matmul(S, cons) f = (DEs[0], DEs[1], DEs[2]) a = jacobian(f, (R, ATP, H) ) g = a.subs(R == x[0], ATP == x[1], H == x[2]) return g
error = 10**(-5)
def backeuler(steps = 100, I=100, Ro = Ro, ATPo = ATPo, Ho = Ho, dt = dt*.01, error = error): x0 = [Ro, ATPo, Ho] a = len(x0) sol = np.zeros((steps, a)) sol[0] = x0 for i in range(steps-1): guess = sol[i] jaco = np.array(dt*dF(guess), dtype = np.float) DE, jpq = Hmodel(sol[i], Ro = Ro, ATPo = ATPo, I=I) #need function dF to evaluate Jacobian at at x xn = guess - np.matmul(inv(np.subtract(eye(a),jaco)),dt*DE) while norm(xn-guess) > error: guess = xn jaco = np.array(dt*dF(guess), dtype = np.float) xdiff = np.subtract(guess, sol[i]) DE, jpq = Hmodel(guess, Ro = Ro, ATPo = ATPo, I=I) xn = guess - np.matmul(inv(np.subtract(eye(a), jaco)),np.subtract(xdiff, dt*DE)) sol[i+1] = xn for k in range(a): y = [x[k] for x in sol] plt.plot(range(steps), y, label="{}".format(vs[k])) print(Hmodel(sol[-1], Ro = Ro, ATPo = ATPo, I=I)) plt.legend() plt.show() backeuler()
#numerical simulation @interact(dt = (.00001, .01, .00001), I=(0,2000,10), steps=(10, 50000, 1), knpq = (.001, 1, .001), kr = (.01, 20000, .01), Ho = (0,100, 1)) def getss(steps = steps, dt = dt, Ro = Ro, Ho=Ho, ATPo=ATPo, PGAo = PGAo, Po=Po, ADPo=ADPo, I=I, O=O, C=C, knpq=.5, kr = kr): x0 = [Ro, ATPo, Ho] a = len(x0) sol = np.zeros((steps, a)) sol[0] = x0 for i in range(steps-1): DE, jpq = Hmodel(sol[i], I=I, knpq=knpq, C=C, O=O, kr =kr, ATPo = ATPo, Ro = Ro, PGAo = PGAo, Po=Po, ADPo=ADPo) step = DE*dt sol[i+1] = np.add(sol[i], step) R = sol[-1][0] fC = kc*C fO = ko*O fR = kr*R jcA = Et/((1/kERC) + (1/(fC+fO)) + 1/fR - 1/(fC+fO+fR)) return jcA
def steadystate(I, C, tsteps=15000, Ro=Ro, ATPo=ATPo, Ho=Ho, PGAo = PGAo, Po=Po, ADPo=ADPo, error = .01, dt = .001): sol= np.zeros((tsteps, 3)) x=sol[0] = [Ro, ATPo, Ho] print(I) for i in range(tsteps-1): DE, jpq = Hmodel(sol[i], I=I, knpq=knpq, C=C, O=O, kr =kr, ATPo = ATPo, Ro = Ro, PGAo = PGAo, Po=Po, ADPo=ADPo) step = DE*dt sol[i+1] = np.add(sol[i], step) R = sol[tsteps-1][0] fC = kc*C fO = ko*O fR = kr*R jcA = Et/((1/kERC) + (1/(fC+fO)) + 1/fR - 1/(fC+fO+fR)) return jcA
@interact() def Isteadystate(Irange = 1000, steps = 10, Ro=Ro, ATPo=ATPo, Ho=Ho, C=C, PGAo = PGAo, Po=Po, ADPo=ADPo, error = .01, dt = .001): Is = np.linspace(.1, Irange, steps) S = [] print('test') for I in Is: jcA = steadystate(I, C, Ro=Ro, ATPo=ATPo, Ho=Ho, PGAo = PGAo, Po=Po, ADPo=ADPo, error = error, dt = dt) S.append(jcA) print(I, jcA) plt.plot(Is, S) plt.show()