Model

Abstract

We used a series of ordinary differential equations (ODEs) to model the ability of the malate sensor to sense different malate concentrations, and a separate series of ODEs to model the ability of the phosphate sensor to sense different phosphate concentrations. The solutions to these equations could predict the luminescence emitted because of this process; the luminescence produced by the malate sensor was predicted to increase for greater initial malate concentrations, whilst the luminescence emitted by the phosphate sensor was predicted to decrease for greater initial phosphate concentrations. The equations were fitted to RLU data and provided a good fit, suggesting our modelling approach was valid. The results showed a greater initial concentration of malate led to larger peaks of RLU, and these peaks occurred earlier in the experiments. These sensors could allow us to infer the activity of PhoBac, as it releases phosphate in the presence of malate and uptakes excess phosphate in waste water.

Introduction

PhoBac has been engineered to release phosphate in the presence of malate, as well as uptake phosphate in wastewater. Two sensors, a malate sensor and a phosphate sensor, have been developed by our iGEM team, that could be incorporated into PhoBac to infer its activity. As it senses a malate concentration, the malate sensor produces a Lux protein that has a luminescence; the greater the amount of malate sensed by the malate sensor, the more Lux protein is produced as a result, and thus the more luminescence is observed. On the other hand, the phosphate sensor would “shut down” in the presence of too much phosphate, and would produce fewer Lux proteins nd, subsequently, emit less luminescence. To better understand these processes and try to predict the performance of the sensors in response to different concentrations of malate and phosphate, we developed two mathematical models, each consisting of a series of ordinary differential equations (hereafter referred to as ODEs). The models give the expected evolution of the amount of luminescence emitted over time; this evolution depends on the initial concentration of malate in the first model, and on the initial concentration of phosphate in the second model. We fit the models to data and subsequently analyse how long it takes for the luminescence emitted by the malate and phosphate sensors to peak, and how large this effect is, and thus how changing the initial concentration of malate or phosphate affects this activity.

Methods

Our modelling approach was adapted from previous models developed by Radeck et al. (2013) [1] and Fritz et al. (2015) [2]. To model the production of Lux protein by the malate sensor in the presence of malate, we considered two ODEs for each initial concentration of malate: one for the concentration of mRNA at time t, m(t), that increases due to transcription and decreases due to decay, and one for the concentration of the Lux protein at time t, p(t), that increases due to translation of the mRNA and decreases due to decay. The units of measurement for translation are set to be RLU/mRNA/min, so that p(t) is changed to be the luminescence released from the malate sensor in the second ODE. Our approach was identical for modelling the activity of the phosphate sensor in the presence of phosphate. We are assuming that only one Lux protein is produced, and that luminescence is directly proportional to the amount produced of this Lux protein. We are also assuming that changes in the luminescence signal are caused by changes in the transcription rate alone. The equations used for each initial concentration of malate were as follows:

Here, α(t) represents the time-dependent apparent transcription rate, τm represents the decay rate of the mRNA, β represents the translation rate and τp represents the decay rate of the Lux protein. The parameters β, τm and τp were taken to be constant throughout the simulation, whilst α(t) was allowed to vary over time. We updated the parameter value for α(t) at every observation; 27 observations were made for each experiment, therefore we used 27 different values for α(t) in the model. Thus, the model effectively had 30 parameters with values to be determined. Our initial estimates were based on those used by Radeck et al. (2013) [1]. As data was unavailable for the mRNA concentrations, the parameter values for the ODEs were chosen to give the best fit to the luminescence data. The ODEs were fit to RLU data, where this data was obtained from observing PhoBac in different initial concentrations of malate over 130 minutes. Luminescence observations were made at the start of the experiment and over 26 subsequent cycles (with each cycle lasting approximately 5 minutes). The experiment was repeated three times for each malate concentration, and the ODEs were fitted to the average RLU observed for each initial malate concentration. Thirteen different initial concentrations of malate were considered, thus we fit the data to thirteen systems of ODEs, with each system comprising the two equations discussed above. The data was fit using the lsqcurvefit function on MATLAB, which solves non-linear curve fitting problems using a least squares approach. Specifically, the Levenberg-Marquardt algorithm was used, as we fit 30 parameters to 27 observations and our system was underdetermined. Levenberg (1944) [3], Marquardt (1963) [4] and Moré (1977) [5] have written in greater detail about this algorithm.

Results

The predicted luminescence outputs for the best-fitting ODEs for each initial concentration of malate are shown in the figure below.

For small initial concentrations of malate (i.e. between 0 and 0.02 nM), the luminescence grows very slowly over time, with the predicted peak being no greater than about 10,000 RLU and occurring at the end of the experiment. As the initial malate concentration grows up to 0.05 nM, the predicted peak luminescence stays approximately the same, but occurs earlier in the experiment (after roughly 120 minutes for 0.03 nM of malate and after roughly 70 minutes for 0.05 nM of malate). At initial concentrations beyond this, the predicted peak luminescence increases, up to a maximum of approximately 60,000 RLU for initial concentrations 1 nM or more. Additionally, the time to reach peak luminescence also initially appears to occur earlier as the initial concentration gets greater, with the earliest predicted time to see peak luminescence being around 45 minutes for an initial malate concentration of 0.3 nM. At initial concentrations greater than this, the predicted time to observe peak luminescence gets progressively later. These results mostly match our expectations, but the fact that the luminescence increases over time for an initial malate concentration of 0 nM is surprising: the malate sensor should not be sensing any malate, and therefore not transcribing mRNA to be translated into Lux protein. However, the luminescence observed is orders of magnitude smaller than that observed for larger initial malate concentrations, so this may simply be due to measurement errors. The predicted luminescence outputs for the best-fitting ODEs for each initial concentration of phosphate are shown in the figure below.

The figures show that, as the initial concentration of phosphate is increased, the maximum luminescence observed is reduced. The phosphate sensor effectively seems to “shut down” in the presence of excessive phosphate.

Conclusion

We have modelled the activity of the malate and phosphate sensors by fitting data from wet lab experiments to a series of ODEs for each set of initial conditions. Our models seem to fit the data reasonably well; our modelling approach thus appears to be valid for this experiment. The results of the first experiment suggest that the malate sensor produces greater amounts of Lux protein in larger initial concentrations of malate, and up to a point, this activity occurs more quickly as the initial malate concentration is increased. Furthermore, the second experiment results show that the phosphate sensor produces less Lux protein and “shuts down” in the presence of greater concentrations of phosphate. This suggests that our sensors could be incorporated into our engineered strain of Basillus subtilis, to monitor its ability to release phosphate in the presence of malate and uptake excess phosphate in waste water.

Appendix: MATLAB Code

The following code was used to generate our model. The data used in this example was one run of the experiment for the malate sensor with an initial malate concentration of 0 nM.

                
                  %% This script takes RLU data and gives the best-fitting parameters for the ODEs model, then plots the predicted versus actual output
%% To try different inputs, edit the Sdata entries
% measurements taken approximately every 5 minutes
Time = [0	4.98576666700000	9.98543333300000	14.9853666700000	19.9853333300000	24.9851166700000	29.9848000000000	34.9848666700000	39.9848333300000	44.9847666700000	49.9845833300000	54.9847500000000	59.9845166700000	64.9844166700000	69.9840500000000	74.9840000000000	79.9840833300000	84.9839333300000	89.9839333300000	94.9837666700000	99.9835500000000	104.983566700000	109.983633300000	114.983250000000	119.983166700000	124.983150000000	129.983050000000];
% observed RLU - this is example data with an initial malate concentration
% of 0 nM
Sdata = [25.33333333	40.66666667	34.66666667	30	33.33333333	40.66666667	32.33333333	51.33333333	59.66666667	52.33333333	68	88	111.3333333	138	135.6666667	234.6666667	252	275	334.6666667	446.3333333	542.6666667	622.3333333	788.3333333	885.6666667	1071.333333	1114.333333	1196.333333];
Sdata = Sdata';
% initial parameter estimates taken from Radeck et al. (2013) 
c0 = [5; 1e04; 4.2; 4.2e-03; 4.2e-03; 4.2e-03; 4.2e-03; 4.2e-03; 4.2e-03; 4.2e-03; 4.2e-03; 4.2e-03; 4.2e-03; 4.2e-03; 4.2e-03; 4.2e-03; 4.2e-03; 4.2e-03; 4.2e-03; 4.2e-03; 4.2e-03; 4.2e-03; 4.2e-03; 4.2e-03; 4.2e-03; 4.2e-03; 4.2e-03; 4.2e-03; 4.2e-03; 4.2e-03];
%% this function fits the ODEs to the data using a least-squares approach, given the initial parameter estimates
lb = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
[pbest,presnorm,presidual,exitflag,output] = lsqcurvefit(@fit,c0,Time,Sdata,lb);
%% this function fits the ODEs to the luminescence data
%% (the mRNA is not observed)
function S = fit(c,t,~)
% initial conditions
% (assumed no mRNA initially)(initial RLU set to average at time t = 0)
y0 = [0;27.87179487];
%% This function solves the ODEs given the initial conditions
[T,Sv] = ode45(@DifEq,t,y0); 
    %% this function gives the ODEs to be solved
    function dS = DifEq(t,y)
        ydot = zeros(2,1);
        % alpha is allowed to vary: it takes a different value for every
        % observation
        alpha = c(4)*(t < 4) + c(5)*(t < 9) + c(6)*(t < 14) + c(7)*(t < 19) + c(8)*(t < 24) + c(9)*(t < 29) + c(10)*(t < 34) + c(11)*(t < 39) + c(12)*(t < 44) + c(13)*(t < 49) + c(14)*(t < 54) + c(15)*(t < 59) + c(16)*(t < 64) + c(17)*(t > 69) + c(18)*(t < 74) + c(19)*(t < 79) + c(20)*(t < 84) + c(21)*(t < 89) + c(22)*(t < 94) + c(23)*(t < 99) + c(24)*(t < 104) + c(25)*(t < 109) + c(26)*(t < 114) + c(27)*(t < 119) + c(28)*(t < 124) + c(29)*(t < 129) + c(30)*(t > 129);
        % m: mRNA concentration
        % p: Lux protein concentration (converted directly to luminescence)
        % mRNA is transcripted and decays
        % Lux protein is translated from mRNA and decays
        % dm/dt = alpha - (log(2)/tau_m) * m
        % dp/dt = beta * m - (log(2)/tau_p) * p
        ydot = [alpha - (log(2)/c(1))*y(1);
                c(2)*y(1) - (log(2)/c(3))*y(2)];
        dS = ydot;
    end
S = Sv(:,2);
end
                
              

References

1. Radeck, Jara, et al. "The Bacillus BioBrick Box: generation and evaluation of essential genetic building blocks for standardized work with Bacillus subtilis." Journal of biological engineering 7.1 (2013): 1-17.
2. Fritz, Georg, et al. "A new way of sensing: need-based activation of antibiotic resistance by a flux-sensing mechanism." MBio 6.4 (2015): e00975-15.
3. Levenberg, Kenneth. "A method for the solution of certain non-linear problems in least squares." Quarterly of applied mathematics 2.2 (1944): 164-168.
4. Marquardt, Donald W. "An algorithm for least-squares estimation of nonlinear parameters." Journal of the society for Industrial and Applied Mathematics 11.2 (1963): 431-441.
5. Moré, Jorge J. "The Levenberg-Marquardt algorithm: implementation and theory." Numerical analysis. Springer, Berlin, Heidelberg, 1978. 105-116.