Team:Lund/Modeling/Implementation

Modeling

Overview

Numerical reliability and reproducibility is a crucial part of research. In this section we provide a detailed description of the implementation of our model.

Programming environment

The entire system was implemented in MATLAB (2017a, The Mathworks Inc.) using the SimBiology© app which provides a block diagram environment for a model-based design implementation of biochemical networks. By using drag and drop features, species and reactions are placed in an editor and connected to yield a model based representation of the system. The representation then generates code which constitutes to the differential equations (or transitions in the case of stochastic modeling) which can be used for simulations. The SimBiology app thus provides an efficient way of building biochemical network models and automatically generate differential equations without the need to manually balance reaction fluxes which otherwise can be a highly error prone process.

Fig. 1 illustrates our model implemented in the graphical user interface. The diagram consists of blue boxes, yellow nodes and arrows interconnecting the components. The boxes represent the interacting species and the nodes the corresponding reactions. Lines pointing from a specie into a node represents a reactant and lines pointing out from a node to a specie represents a product. Dashed lines indicate that the reactant is not consumed under the reaction.

The illustrated system can be subdivided into five parts. The first part (a) illustrates the lacI dynamics which expresses the lac repressor. In (b), the IPTG induction of the repressed genes represented by their operons ON and OE is shown. These operons regulate the epression of GFP1-9 and GFP10-ER-GFP11 respectively. The NahR pathway is marked by (c) and shows the binding to NahR to the Psal promoter followed by the subsequent transcriptional activation by salicylate to yield GFP1-9. In (d), GFP10-ER-GFP11 is translated and activates by the addition of estrogen. In the final part (e), the GFP association occurs, yielding a fully restored fluorophore.

Figure 1: Illustration of the SimBiology© implementation of our system. The networks can be considered as five parts (a)-(e) governing the dynamics of separate components. Those include the expression of the lac repressor, lac regulation, NahR and GFP10-ER-GFP11 pathway and final GFP association respectively.

Experimental setup

Simulations

To properly emulate the working conditions of the cell at the time of addition in the wet lab, the system was run until steady state without the addition of inducers. The steady state values were then saved and a new simulation was started using the steady state values as initial concentrations and adding IPTG, salicylate and estrogen to induce gene expression. The inital simulation was run with initial concentrations set to zero expect for the gene content which was set at 83.2 nM for $\ce{O_N}$, $\ce{O_E}$ and $\ce{I_R}$ (see parameter section for details).

Calculation of rise time and settling time

To measure the responsiveness of the biosensor we used the rise time and settling time as measures. The rise time is defined as the time it takes for a signal to go from 10% to 90% of the steady state value. The settling time is defined at the time it takes for a signal to fall within 2% of the steady state value. These entities were calculated using the function STEPINFO (MATLAB 2017a, The Mathworks Inc.).

Calculation of global sensitivities

The global sensitivities were calculated using the Sensitivity Analysis Library (SALib, see [1]) using Python 3.5.3. The method used was the Sobol method which is a variance based sensitivity analysis method. This was chosen because it has been recommended to use variance based methods whenever the computational cost is not a problem [2]. The input parameters were chosen by logarithmically sampling the parameter space using Saltelli's sampling scheme. For more details, see the documentation page.

Implementation details

Choice of ODE solver

Since we had parameters differencing with several orders of magnitude, numerical stiffness was expected and a stiff ordinary differential equation (ODE) solver had to be used. For this we choose MATLABs solver ODE23t (MATLAB 2017a, The Mathworks Inc.) as it was found to have the best performance on our system among the other solvers in the MATLAB ecosystem. This solver is commonly used to solve moderately stiff ODEs and uses a variant of the classical trapezoidal rule.

Choice of optimization algorithm

In performing the parameter estimation, several optimization algorithms were evaluated and tried such as FMINSEARCH (MATLAB 2017a, The Mathworks Inc.), LSQNONLIN (MATLAB 2017a, The Mathworks Inc.) and PARTICLESWARM (MATLAB 2017a, The Mathworks Inc.). The function LSQNONLIN was chosen as it seemed to perform best. The solver is commonly used to solve nonlinear least square problems and uses the Levenberg-Marquart algorithm. All optimizations were run using standard configurations.

Simulation parameters

To ensure numerically reliable solutions, all simulations were carried out using an absolute error tolerance of 10-8 and carefully inspected to keep off spurious solutions. In addition, several solutions were tested against other ODE solvers such as ODE15s (MATLAB 2017a, The Mathworks Inc.) to ensure reliability.

References

  1. [1] Herman, J. and Usher, W. (2017) SALib: An open-source Python library for sensitivity analysis. Journal of Open Source Software, 2(9).
  2. [2] Zi, Z. (2011). Sensitivity analysis approaches applied to systems biology models. IET Systems Biology, 5(6), pp.336-346.