In this Animation we can see a Schematic Visualistion of the Process.
First, the wastewater that transports phosphate is filled in the reactor with flowrate F1, or flows throuh the reactor continuely. While the light is on, our PBP binds the phosphate and therefore seperates it from the wastewater particles. Since they are retained by a membrane, the PBPs cannot leave the reactor. It can happen that wastewater particles stuck at the membrane, but this can be avoided by cleaning the reactor occasionally.
When most of the PBPs have bounded some phosphate, the inflow is stopped and the remaining wastewater leaves the reactor. Now we fill it with a solution and switch off the light inside the reactor. The phosphate is then released into the solution and can be extracted with flow rate F2. From this solution you can extract the phosphate and process it to polyphosphate.
The light in the reactor is than activated again and after some time the PBPs can be used again.
Here $k_1$ indicates the rate of how fast phosphate is bounded by the modified PBP. $k_2$ gives the information about the change from active PBP to inactive PBP induced by light and k3 stands for the phosphate release rate. The regeneration of the inactive PBP back to the active state is shown by $k_4$. Finally the back reactions $k_{-1},k_{-2}$ and $k_{-4}$ depend on the strength of the forward reactions. The back reaction of $k_3$ is not needed, since inactive PBP is not able to bind phosphate.
The black and orange arrows indicate reactions that happen in light or in darkness, and are later modeled by a binary switch.
First some Notation:
By using Michaelis-Menten-Kinetic we can write the reactions in a more compact form, where $v_1,\dots ,v_{-4}$ denote the speed of the reaction.
Now the Speed $v_i$ can be specified with the rate of change of concentration ($k_i$).
A reaction $S_1 + S_2 \overset{k_1}{\rightarrow} P$ is written as:
$$ \frac{\partial P}{\partial t} = v = k_1\cdot S_1 \cdot S_2 $$In conclusion we get:
$$\frac{\partial}{\partial t} \begin{pmatrix} PhA \\ Ph \\ Ep \\ En \\ Epp \\Enp \end{pmatrix}= \begin{pmatrix} -L*k_1*Ep*PhA & L*k_{-1}*Epp & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & k_3*Enp*(1-L) & 0 & 0\\ -L*k_1*PhA*Ep & L*k_{-1}*Epp & 0& 0& 0& L*k_4*En & -(1-L)*k_{-4}*Ep \\ 0& 0& 0& 0& k_3*Enp*(1-L)& -L*k_4*En & (1-L)*k_{-4}*Ep \\ L*k_1*Ep*PhA & -L*k_{-1}*Epp& -(1-L)*k_2*Epp& L*k_{-2}*Enp &0 &0 &0 \\ 0&0 & (1-L)*k_2*Epp & -L*k_{-2}*Enp &-k_3*Enp*(1-L) & 0 &0 \end{pmatrix} $$Therefore we obtain a system of ordinary differential equations.
$$ \partial PhA=-L*k1*Ep*PhA+L*k1n*Epp $$$$ \partial Ph=k3*Enp*(1-L) $$$$ \partial Ep=-L*k1*Ep*PhA+L*k1n*Epp+k4*En*L-k4n*Ep*(1-L) $$$$ \partial En=k3*Enp*(1-L)-k4*En*L+k4n*Ep*(1-L) $$$$ \partial Epp=k1*Ep*PhA*L-k2*Epp*(1-L)+k2n*Enp*L-k1n*Epp*L $$$$ \partial Enp=k2*Epp*(1-L)-k2n*Enp*L-k3*Enp*(1-L) $$Next, we want to include the flows. Since only phosphate wastewater and phosphate are affected by the inflows and outflow, these are the only ones we have to adapt. The mass bilance is determined by the change of concentration. Since it is difficult to add up concentrations, you can formulate it with flows: $$(C\cdot V)'=V\cdot C' = f_{In}\cdot C_{In}- f_{Out}\cdot C$$ With the assumption that the inflow $f_{In}$ equals the outflow $f_{Out}$ and is denoted by $f$, you get the following formula: $$C'=\frac{f}{V}\cdot (C_{In}-C)$$ Here $C_{In} $ denotes the starting concentration of the product, that get filled in the reactor.
$F1$ is now the flow of the wastewater through the reactor and $F2$ is the flow of the solution, that transports the released Phosphate after the Light is off.
We obtain the following System of ODEs:
$$ \partial PhA=-L*k1*Ep*PhA+L*k1n*Epp + \frac{F1}{V}(PhA_{In}-PhA)*L$$$$ \partial Ph=k3*Enp*(1-L) - \frac{F2}{V}Ph*(1-L) $$$$ \partial Ep=-L*k1*Ep*PhA+L*k1n*Epp+k4*En*L-k4n*Ep*(1-L) $$$$ \partial En=k3*Enp*(1-L)-k4*En*L+k4n*Ep*(1-L) $$$$ \partial Epp=k1*Ep*PhA*L-k2*Epp*(1-L)+k2n*Enp*L-k1n*Epp*L $$$$ \partial Enp=k2*Epp*(1-L)-k2n*Enp*L-k3*Enp*(1-L) $$For simplicity we set $V=1$ since this is just the volume scale of the reactor and does not influence the reaction-scheme.
After that we implement the modell as the following:
import numpy as np
def Kinetic(t,z,k1,k1n,k2,k2n,k3,k4,k4n,L,tOff,tOn,F1=0,F2=0,PhA_In=0):
PhA,Ph,Ep,En,Epp,Enp = z
dPhA=-k1*Ep*PhA*L(t,tOff,tOn)+k1n*Epp*L(t,tOff,tOn) + F1*(PhA_In-PhA)*L(t,tOff,tOn)
dPh=k3*Enp*(1-L(t,tOff,tOn)) - F2*Ph*(1-L(t,tOff,tOn))
dEp=-k1*Ep*PhA*L(t,tOff,tOn)+k1n*Epp*L(t,tOff,tOn)+k4*En*L(t,tOff,tOn)-k4n*Ep*(1-L(t,tOff,tOn))
dEn=k3*Enp*(1-L(t,tOff,tOn))-k4*En*L(t,tOff,tOn)+k4n*Ep*(1-L(t,tOff,tOn))
dEpp=k1*Ep*PhA*L(t,tOff,tOn)-k2*Epp*(1-L(t,tOff,tOn))+k2n*Enp*L(t,tOff,tOn)-k1n*Epp*L(t,tOff,tOn)
dEnp=k2*Epp*(1-L(t,tOff,tOn))-k2n*Enp*L(t,tOff,tOn)-k3*Enp*(1-L(t,tOff,tOn))
return np.array([dPhA,dPh,dEp ,dEn ,dEpp ,dEnp ])
L is here the light function. It is 1, when the light is on and 0 if the light is off.
$$L(t):= \begin{cases} 1 & \text{if Light on} \\ 0 & \text{if Light off} \end{cases} $$#Light On and Off:
def L(t,tOff,tOn):
if t < tOff:
return 1
else:
if t< tOn:
return 0
else:
return 1
For the solution we used an implizit ODE-Solver, since our system ist stiff. The stiffness is shown in A1. Explicit we used the scipy.integrate solver with the Radau-method, which is a implicit Runge-Kutta method of Radau IIA family of order 5.
We switch off the light at the time tOff given in minutes.
from scipy.integrate import solve_ivp
def Solution(k1,k1n,k2,k2n,k3,k4,k4n,L,tOff,tOn,length,start,F1=0,F2=0):
#k1,k1n,k2,k2n,k3,k4,k4n,L, t for Light Off, t for Light On, F1, F2, Start
parameter=(k1,k1n,k2,k2n,k3,k4,k4n,L,tOff,tOn,F1,F2,start[0])
#Solution
sol= solve_ivp(Kinetic,(0,length),start,method="Radau", args=parameter)
return sol.y, sol.t
Lets look at the solution of this modell for some parameters. The parameters are orientated at the process times of the enzymes, that we used (further informations in A2). But they can change by using diffrent PBP or lightswitches. The concentration of each component is given in g/l. We can exchange the parameters easily and look at the change in the solution.
#Parameter
k1=0.024
k1n=0.002
k2=0.5
k2n=0.2
k3=0.024
k4=0.0017
k4n=0.0084
tOff=30 # minutes
#Initial Conditions:
#PhA,Ph,Ep,En,Epp,Enp
start=[100,0,100,0,0,0]
length=300
from ODE_Plots_iGEM import Plot_ODE_All
y,t=Solution(k1,k1n,k2,k2n,k3,k4,k4n,L,tOff,length,length,start)
Plot_ODE_All(y,t)
Now we want to include the flows and take a look what happens for diffrent flows and tOFFs.
from ODE_Plots_iGEM import Subplot_4
start_Flow=[10,0,100,0,0,0]
y1,t1=Solution(k1,k1n,k2,k2n,k3,k4,k4n,L,tOff,length,length,start_Flow,2) #Big Inflow
y2,t2=Solution(k1,k1n,k2,k2n,k3,k4,k4n,L,tOff,length,length,start_Flow,0.01) #Small Inflow
y3,t3=Solution(k1,k1n,k2,k2n,k3,k4,k4n,L,10,length,length,start_Flow,0.2) #Short Light On - 10 Minutes
y4,t4=Solution(k1,k1n,k2,k2n,k3,k4,k4n,L,200,length,length,start_Flow,0.2) #Long Light On - 200 Minutes
Subplot_4(y1,y2,y3,y4,t1,t2,t3,t4)
We see that when the flow is high, a lot of phosphate is recovered, but a lot of phosphate is also lost in the process. However, if the flow is too low, not enough phosphate is bound, and you only get a small yield.
The time at which the light is turned off and directed into the harvesting process is also important. If this time is too short, not enough phosphate is bound and the yield is small. If you wait too long, all the enzymes are utilized and you waste unnecessary time.
The time while the light is turned off and directed into the harvesting process is also important. If this time is too short, not enough phosphate is bound and the yield is small. If you wait too long, all the enzymes will be used up and you will lose time unnecessarily.
It is necessary to find the optimal flow and time.
First, we defined an objective function that returns the phosphate production at the end of the interval for a given flow at the optimal time to turn off the light.
The best switch-off time is determined by a program that you can find in the Jupyter Live Script. Here we return the time where the bound active PBP does not change significantly anymore. Therefore, $$\frac{\partial Epp}{\partial t} < TOL \approx 0$$ So the change is smaller than a tolerance TOL which is choosen close to zero.
def Best_Time_at_Flow(F,k1,k1n,k2,k2n,k3,k4,k4n,L,length,start,TOL=0.01):
#PA,P,Ep,En,Epp,Enp
sol_y,sol_t=Solution(k1,k1n,k2,k2n,k3,k4,k4n,L,length,length,length,start,F,0)
index=-1
Epp_Max=-1
for i in range(np.size(sol_y[0])):
l=L(sol_t[i],length,length)
dEPP=k1*sol_y[2][i]*sol_y[0][i]*l-k2*sol_y[4][i]*(1-l)+k2n*sol_y[5][i]*l-k1n*sol_y[4][i]
if abs(dEPP) < TOL:
index=sol_t[i]
Epp_Max=sol_y[4][i]
return index, Epp_Max # Best Time, Epp at this Time
Now we can define the objective Function. We negate the solution, since we will later use a minimizing optimization, but need the maximum.
def objective(F,k1,k1n,k2,k2n,k3,k4,k4n,L,length,start):
i,maxEpp=Best_Time_at_Flow(F,k1,k1n,k2,k2n,k3,k4,k4n,L,length,start)
sol_y,sol_t=Solution(k1,k1n,k2,k2n,k3,k4,k4n,L,i,length,length,start,F,0)
maxPh=(-1)*sol_y[1][-1]
return maxPh
Now we can use Nelder-mead optimization or a diffrent one to gain the optimum flow and the best time to shut down the light at this flow.
from scipy.optimize import minimize
from scipy.optimize import Bounds
from numpy.random import rand
def Optimum_F1_tOff(r_min,r_max,k1,k1n,k2,k2n,k3,k4,k4n,L,length,start,output=True):
# define range for input
boundsFlow=Bounds(r_min,r_max)
# define the starting point as a random sample from the domain
pt = r_min + rand(1) * (r_max - r_min)
# perform the search
result = minimize(objective, pt, args=(k1,k1n,k2,k2n,k3,k4,k4n,L,length,start), method='nelder-mead',bounds=boundsFlow)
# evaluate solution
BestFlow = result['x']
BestTime, Max = Best_Time_at_Flow(BestFlow,k1,k1n,k2,k2n,k3,k4,k4n,L,length,start)
# summarize the result
if output==True:
print('Status : %s' % result['message'])
print('Total Evaluations: %d' % result['nfev'])
print('Best Flow:', BestFlow)
print('Best Time for Light Off:', BestTime)
return BestTime, BestFlow
Since the results may vary, due to the stochastric nature of the algorithm, we let it run a few times and take the mean value.
def Mean_Run(r_min,r_max,k1,k1n,k2,k2n,k3,k4,k4n,L,length,start,n_iter=10):
Flow=[]
for i in range(n_iter):
T,F=Optimum_F1_tOff(r_min,r_max,k1,k1n,k2,k2n,k3,k4,k4n,L,length,start,False)
Flow.append(F[0])
Flow_mean=np.mean(Flow)
Time_mean,Max=Best_Time_at_Flow(Flow_mean,k1,k1n,k2,k2n,k3,k4,k4n,L,length,start)
print('Mean of Best Flows:', Flow_mean)
print('Best Time for Light Off at this Flow:', Time_mean)
return Time_mean, Flow_mean
#Single Run
#BestTime, BestFlow = Optimum_F1_tOff(0,0.9,k1,k1n,k2,k2n,k3,k4,k4n,L,length,start_Flow)
BestTime,BestFlow = Mean_Run(0,0.9,k1,k1n,k2,k2n,k3,k4,k4n,L,length,start_Flow,5)
#Return Solution at the Optimum
y_opt,t_opt=Solution(k1,k1n,k2,k2n,k3,k4,k4n,L,BestTime,length,length,start_Flow,BestFlow,0)
Plot_ODE_All(y_opt,t_opt)
C:\Users\Asus 8600k\AppData\Local\Temp\ipykernel_8944\924294074.py:13: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray. return np.array([dPhA,dPh,dEp ,dEn ,dEpp ,dEnp ])
Mean of Best Flows: 0.8153654551344595 Best Time for Light Off at this Flow: 48.016883762284124
Finally lets have a look what happens after we put the light back on:
y_all,t_all=Solution(k1,k1n,k2,k2n,k3,k4,k4n,L,BestTime,200,500,start,0,0)
Plot_ODE_All(y_all,t_all)
When the light is back on, it will take some time to reactivate the PBP. This depends on how fast the regeneration time of the used PBP with the optogenetic switch is.
This could be a problem in the industry. Therefore, we provide a solution to avoid this problem.
The wastewater flows inside reactor 1, where phosphate binds under the influence of light. After about an hour, the flow is directed to the reactor 2 and in reactor 1 flows solution. While in R2 phosphate is bound, phosphate is released in R1 in the dark state and therefore provides usable phosphate. After the next hour, R1 is emptied and the light is turned on again. Now the PBPs can return to the active state, while in R2 phosphate is released and in R3 phosphate is bound. After the regeneration time in R1, the fusion protein can be used again and we obtain a cycle that gives us a maximum amount of phosphate.
import SALib
from SALib.sample import saltelli
from SALib.analyze import sobol
from scipy.integrate.odepack import odeint
def Kinetic_Sobol(z,t,k1,k1n,k2,k2n,k3,k4,k4n,L,tOff,tOn,F1=0,F2=0,PhA_In=0):
PhA,Ph,Ep,En,Epp,Enp = z
dPhA=-k1*Ep*PhA*L(t,tOff,tOn)+k1n*Epp*L(t,tOff,tOn) + F1*(PhA_In-PhA)*L(t,tOff,tOn)
dPh=k3*Enp*(1-L(t,tOff,tOn)) - F2*Ph*(1-L(t,tOff,tOn))
dEp=-k1*Ep*PhA+k1n*Epp+k4*En*L(t,tOff,tOn)-k4n*Ep*(1-L(t,tOff,tOn))
dEn=k3*Enp*(1-L(t,tOff,tOn))-k4*En*L(t,tOff,tOn)+k4n*Ep*(1-L(t,tOff,tOn))
dEpp=k1*Ep*PhA*L(t,tOff,tOn)-k2*Epp*(1-L(t,tOff,tOn))+k2n*Enp*L(t,tOff,tOn)-k1n*Epp
dEnp=k2*Epp*(1-L(t,tOff,tOn))-k2n*Enp*L(t,tOff,tOn)-k3*Enp*(1-L(t,tOff,tOn))
return np.array([dPhA,dPh,dEp ,dEn ,dEpp ,dEnp ])
def Sobol_analysis():
length=200
y_0=[10,0,100,0,0,0]
t=np.linspace(0,200,400)
#y_sens,t_sens=Solution(k1,k1n,k2,k2n,k3,k4,k4n,L,60,length,length,y_0,0.5,0)
Solution_y=odeint(Kinetic_Sobol,y_0, t, args=(k1,k1n,k2,k2n,k3,k4,k4n,L,60,length,0.5,0))
problem={
'num_vars': 13,
'names': ['k1','k1n','k2','k2n','k3','k4','k4n','PhA_0','Ph_0','Ep_0','En_0','Epp_0','Enp_0'],
'bounds': [[0,0.1], [0,0.1],[0,1], [0,1], [0,0.1],[0,0.1],[0,0.1],[5,20], [0,1],[90,120],[0,1],[0,1],[0,1]]
}
# Generate samples
vals = saltelli.sample(problem, 100)
# initializing matrix to store output
Y = np.zeros([len(vals),1])
# Run model (example)
# numerically soves the ODE
# output is X1, X2, and X3 at the end time step
Y = np.zeros([len(vals),6])
for i in range(len(vals)):
Y[i][:] = odeint(Kinetic_Sobol,[vals[i][7],vals[i][8],vals[i][9],vals[i][10],vals[i][11],vals[i][12]],t, args=(vals[i][0],vals[i][1],vals[i][2],vals[i][3],vals[i][4],vals[i][5],vals[i][6],L,60,length,length,0.5,0))[len(Solution_y)-1]
# completing soboal analysis for each X1, X2, X3, X4,X5,X6
print('\n\n====PhA Sobol output====\n\n')
Si_X1 = sobol.analyze(problem, Y[:,0], print_to_console=True)
print('\n\n====Ph Sobol output====\n\n')
Si_X2 = sobol.analyze(problem, Y[:,1], print_to_console=True)
print('\n\n====Ep Sobol output====\n\n')
Si_X3 = sobol.analyze(problem, Y[:,2], print_to_console=True)
print('\n\n====En Sobol output====\n\n')
Si_X4 = sobol.analyze(problem, Y[:,3], print_to_console=True)
print('\n\n====Epp Sobol output====\n\n')
Si_X5 = sobol.analyze(problem, Y[:,4], print_to_console=True)
print('\n\n====Enp Sobol output====\n\n')
Si_X6 = sobol.analyze(problem, Y[:,5], print_to_console=True)
#Sobol_analysis() #takes a while
In the tables one can read the sensitivity of the individual parameters. Here, S1 stands for the 1st order sensitivity and S1conf for the corresponding confidence interval. ST stands for the total order sensitivity and STconf for the corresponding confidence interval. The 2nd order sensitivity indicates how the output behaves when two parameters interact and the total order considers the effect of all parameter interactions.
We can see that the Parameter $k_{-1}$ is the most sesnsitive for PhA, Ph,Epp and Enp and $k_{-4}$ is the most sensitive for Ep and En.
Now we can determine the multidimensional derivative of $f$ for specific point $x_0$: $$Df(y_0) = (\frac{\partial f(y_0)}{\partial PhA},\frac{\partial f(y_0)}{\partial Ph},\frac{\partial f(y_0)}{\partial Ep},\frac{\partial f(y_0)}{\partial En},\frac{\partial f(y_0)}{\partial Epp},\frac{\partial f(y_0)}{\partial Enp})$$
$$=\begin{pmatrix} -k1*Ep*L & 0 & -k1*PhA*L & 0 & k1n*L & 0 \\ 0 & 0 & 0 & 0 & 0 & k3 \\ -k1*Ep*L & 0 & -k1*PhA*L-k4n*(1-L) & k4*L & k1n*L & 0 \\ 0 & 0 & k4*(1-L) & -k4*L & 0 & k3 \\ k1*Ep*L & 0 & k1*PhA*L & 0 & -k2*(1-L)-k1n*L & k2n*L \\ 0 & 0 & 0 & 0 & k2*(1-L) & -k2n*L-k3 \end{pmatrix}$$Therefore $f(y)$ ist continous and differentiable. Since the individual parameters $k1,k_{-1}, k2,k_{-2},k_3,k_4,k_{-4}$ are fixed, as well as $y_0$ lie only on a restricted range, the derivation is additionally restricted. Thus the global Lipschitz continuity [https://en.wikipedia.org/wiki/Lipschitz_continuity] of the function follows. By the theorem of Picard-Lindelof [https://en.wikipedia.org/wiki/Picard%E2%80%93Lindel%C3%B6f_theorem] it holds that locally a unique solution of the differential equation exists which is continuous.
A initial value problem $y'=f(t,y)$ with $y(t_0)=y_0$ is stiff for a solution $y(t)$, if the eigenvalues $\lambda(t)$ of the jacobian matrix $Df$ fullfill $$ \kappa (t):=\frac{\max_{Re(\lambda (t))<0}|Re \lambda(t)|}{\min_{Re(\lambda (t))<0}|Re \lambda(t)|} >> 1$$ Here $\kappa (t)$ is the stiffness rate.
Stiff problems have big diffrences in their decay behavior for each component. A explicit solver would need to much steps to solve such a problem.
Therefore, if a problem is stiff it is highly recomended to use an implicit solver.
We will test the stiffness for the choosen parameters and initial vector, with $F1=F2=0$ and $tOff=30$.
from ODE_Plots_iGEM import array_vektor
lenght=200
par=[k1,k1n,k2,k2n,k3,k4,k4n,1]
y_0=[100,0,100,0,0,0]
y_stiff,t_stiff=Solution(k1,k1n,k2,k2n,k3,k4,k4n,L,30,length,length,y_0,0,0)
y_vec = array_vektor(y_stiff)
def Eigenvalues(par, p,t):
l=L(t,tOff,length)
PhA=p[0]
Ph=p[1]
Ep=p[2]
En=p[3]
Epp=p[4]
Enp=p[5]
Jacobi=np.array([[-k1*Ep*l,0,-k1*PhA*l,0,k1n*l,0],[0,0,0,0,0,k3],[-k1*Ep*l,0,-k1*PhA*l-k4n*(1-l),k4*l,k1n*l,0],[0,0,k4*(1-l),-k4*l,0,k3], [k1*Ep*l,0,k1*PhA*l,0,-k2*(1-l)-k1n*l,k2n*l],[0,0,0,0,k2*(1-l),-k2n*l-k3]])
w, v = np.linalg.eig(Jacobi)
return w
def stiff(eigenvalue):
EW_neg=[]
for i in eigenvalue:
real=i.real
if real<0:
EW_neg.append(abs(real))
EW_max=np.max(EW_neg)
EW_min=np.min(EW_neg)
stiff_fac=EW_max/EW_min
return stiff_fac
def Run():
List=[]
index=0
for i in range(len(y_vec)):
t=t_stiff[i]
y=y_vec[i]
value= Eigenvalues(par, y,t)
value_stiff=stiff(value)
List.append(value_stiff)
#print([index,value,value_stiff]) #print for the big Output
#print(value_stiff, sep = ", ") #unsorted
index+=0.2
List.sort()
print(List)
return List
Run()
[59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 131.76470588235296, 131.76470588235296, 131.76470588235296, 131.76470588235296, 131.76470588235296, 131.76470588235296, 131.76470588235296, 131.76470588235296, 138.64493492474207, 420.4042336256544, 676.4037106448961, 1100.060821943348, 1723.3729764088894, 1.9374693750321764e+16, 2.2153351674846536e+16, 2.2834053700227012e+16, 2.723495081958813e+16, 7.667614907898059e+16, 7.750968630503605e+16, 1.437574435673917e+17, 1.4860299054609754e+17, 2.094759212088524e+17, 2.681972413147013e+17, 3.07116562596575e+17, 3.9013556145594726e+17, 4.857179050041234e+17, 5.081205684601648e+17, 7.32687994379471e+17, 7.979805995955448e+17, 1.0690280115952003e+18, 1.142988864305846e+18, 1.4093917784950917e+18, 1.8276056331952952e+18, 1.8689349658458168e+18, 2.999878887988493e+18, 3.7518839440333906e+18, 3.983880969189397e+18, 4.1611350263806e+18, 5.90625335791317e+18, 5.993959261187721e+18, 6.971284736905508e+18, 7.540310095143675e+18, 8.224899983603799e+18, 1.3780182548198556e+19, 1.4580634449677543e+19, 2.4524443173870453e+19]
[59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 59.523809523809526, 131.76470588235296, 131.76470588235296, 131.76470588235296, 131.76470588235296, 131.76470588235296, 131.76470588235296, 131.76470588235296, 131.76470588235296, 138.64493492474207, 420.4042336256544, 676.4037106448961, 1100.060821943348, 1723.3729764088894, 1.9374693750321764e+16, 2.2153351674846536e+16, 2.2834053700227012e+16, 2.723495081958813e+16, 7.667614907898059e+16, 7.750968630503605e+16, 1.437574435673917e+17, 1.4860299054609754e+17, 2.094759212088524e+17, 2.681972413147013e+17, 3.07116562596575e+17, 3.9013556145594726e+17, 4.857179050041234e+17, 5.081205684601648e+17, 7.32687994379471e+17, 7.979805995955448e+17, 1.0690280115952003e+18, 1.142988864305846e+18, 1.4093917784950917e+18, 1.8276056331952952e+18, 1.8689349658458168e+18, 2.999878887988493e+18, 3.7518839440333906e+18, 3.983880969189397e+18, 4.1611350263806e+18, 5.90625335791317e+18, 5.993959261187721e+18, 6.971284736905508e+18, 7.540310095143675e+18, 8.224899983603799e+18, 1.3780182548198556e+19, 1.4580634449677543e+19, 2.4524443173870453e+19]
We see that the smallest stiffness rate is about 19 and than increases highly. Therefore we have indeed a stiff problem and need an implicit solver.
For a first-order differential equation of a decreasing behavior $$u(t)'=-k\cdot u(t) $$ with initial value $u(t_0)=u_0 we get the solution: $$u(t)=u_0\cdot e^{-k\cdot t} $$
Now we can solve this for the unknown parameter $k$. $$ k=\frac{\log(\frac{u(t)}{u_0})}{t} $$
Lets take a look at one of the reactions: $$Ep + PhA \leftrightarrow Epp $$
Lets say the half-life time is at approx. $t=12$ for this type of enzyme.
Than the equation reduces to the following: $$ k=\frac{\log(\frac{1}{2})}{15} = -0.025 $$
Therefore $k_1\approx 0.025$
The backward reaction $k_{-1}$ can be determined by the steady state and strength of binding.
In the steady state we have $Ep'=0$ and conclusive $$k_1\cdot Ep \cdot PhA = k_{-1}Epp .$$ Or $$k_{-1}=k_1\cdot \frac{Ep\cdot PhA}{Epp} $$
If the phosphate is strongly bounded to Ep than the quotient $\frac{Ep\cdot PhA}{Epp}$ should be small (below one). If the binding is weak the quotient is big.
In our case we have a quite strong binding. So we obtained for $k_{-1}=0.002$
This values will be by far not accurate, but at the moment this is the best we can do for approximating them. You could improve this if you have good literatur values or detailed lab reports.