Model Description

In this model, we included transcription, translation, translocation through E. coli’s membrane, protein folding and mRNA and protein degradation in the cytoplasm and medium. NGF synthesis is put under a T7 promoter repressed by LacI, so we also modelled the IPTG induction, an analogue of lactose (allolactose). Finally, NGF is secreted to the medium through a Type I secretion system in which the export signal peptide is not cleaved during translocation. Our Biobrick is design to synthetize and export TEV protease in order to cleave the signal peptide and thus produce functional NGF.

The molecular mechanisms included in our model appear in the following scheme:

Our model includes the following variables:

The synthesis of NGF and TEV is placed under the control of the Plac (T7) promoter. The promoter can be in two different states: off (Po) when the lacO site is occupied by the repressor lacI, preventing RNA polymerase from binding and thus preventing transcription, or free (P_f) thanks to IPTG binding to the repressor, and thus freeing the operator site. We assume that one IPTG molecule binds with one repressor molecule (one dimer per site, with two sites dimerizing, therefore 4 lactose or IPTG molecules required), freeing the promoter and restoring RNA polymerase binding capacity.

The real mechanism of promoter Plac is more complex, as described in ^[1], but this simplification is sufficient for our model.

The transport of IPTG from outside the cell to cytoplasm is considered to be only due to free diffusion through the membrane by two first order reaction with the same kinetic constant. IPTG is not considered to be degraded neither in the cytoplasm nor in the medium.

For the TEV and NGF transcription, we use a first-order reaction where the rate of mRNA production (m) depends on the concentration of the free promoter (P_f).

For the TEV and NGF translation, we first consider binding of ribosomes to ribosome binding site (the same association constant is used since the RBSs are the same), and then translation rate is proportional to the protein length. Since TEV and NGF have approximately the length, we consider only one translation rate β.

Even though it still has an export peptide, TEV is assumed to be functional in the cytoplasm (although less functional than if it had no export peptide). Since NGF has TEV cleaving site between the coding sequence and the export peptide, a fraction of NGF is cleaved inside the cytoplasm and thus cannot be secreted.

We use a simple model to simulate TEV kinetics: TEV recognizes the signal sequence ENLYFQ, bind to its substrate and then cleaves the export peptide. This process can thus be modeled by the following equations:

K₁, k_-1 and k₂ are taken lower than constants found in literature, in order to model the fact that TEV still has its signal peptide and is consequently less functional than usually.

The transport of NGF and TEV with their export signal peptide from inside the cell to the medium is assumed to follow Michaelis-Menten enzymatic kinetics in which the transporter channel (composed of HlyB in the inner membrane, bound to HlyD and recruiting TolC in the outer membrane) plays the role of the enzyme and intracellular protein the role of the substrate.

Each protein (NGF and TEV) via its export signal peptide HlyA can bind to the HlyB-HlyD complex pore, forming a protein-transporter complex (NGFt or TEVt). Translocation correspond to the dissociation of this complex, resulting in restoring a free transporter and secreting NGF or TEV in the medium (NGF_um and TEV_m), which stand for the products.

This model is valid for one bacterial cell, but for our model to fit with our proof of concept system, which is a microfluidic chip chamber containing 100 μL of bacterial culture, we need to integrate the number of bacteria contained in the chamber. Therefore, our model helps to determine which is the most accurate bacteria amount required to produce the appropriate NGF concentration.

Since in the lab, our easiest way to determine number of bacteria is measuring OD600, we determined which is the best OD600 for the final proof-of-concept experiment, considering that an OD600 value of 1 corresponds to 8.10⁸ bacteria units.

Once in the medium, both NGF and TEV are still bound to the export signal peptide HlyA. We assume there is a very small amount of functional TEV, that is sufficient to cleave TEV signal peptide, producing more functional TEV.

As for the transporter, we use a simple model in which TEV recognizes the signal sequence ENLYFQ, bound to its substrate (which can be either NGF with its export peptide or another TEV molecule with its own export peptide) and then cleaves the export peptide. This process can thus be modeled by the following equations:

Finally, in the cytoplasm and in the medium, mRNAs and proteins are degraded and all degradations are assumed to follow first-order kinetic reactions.

MODEL PARAMETRISATION

From these equations, we obtained a system of differential equations mostly based on mass action kinetics (Ref.)(GET IT HERE). We numerically solved the ordinary differential equations system using Euler method implemented in Python. The constants we used were mainly determined from literature AND are given in table:

We determined the temporal evolution of secreted NGF concentration in the medium, in order to get the u(0,t) term used in our following diffusion model.

The aim of this first model is to demonstrate that we can expect an appropriate secreted recombinant NGF concentration to observe neurite growth. However, we had to make several assumptions to parametrize the model. We scanned different parameter values for the values we assumed (such as number of transporters or kinetic parameters for translocation) in order to check the range of NGF amount we can reasonably expect. We also studied influence of IPTG induction and number of bacteria, since they are parameters our wetlab can control to best fit recombinant NGF secretion with what we need.

We co-transformed our bacteria with a plasmid expressing HlyB and HlyD, two of the components of the secretion pore. However, we did not quantify the number of pores each cell contains, and we are only able to estimate it, based on assumptions made in^[5]. Consequently, we scanned a range of different values for the quantity of transporters in order to see the range of NGF concentration we can expect.

The following graph shows the predicted NGF concentration in the microfluidic chip chamber for a number of pores varying: no pore (A.), 10 per cell (B.), 100 per cell (C.) and 500 per cell (D.):

k₄ estimated value does have a real impact on the total amount of secret NGF, as evidenced by the graph (given here for an IPTG induction level of 1 mM, and 500 transporters/cell). Indeed, when k4 varies from 0.5 min^-1 to 4 min^-1, the amount of total secreted NGF doubles from 60 ng/mL to 120 ng/mL.

As expected the final NGF concentration (both in the cytoplasm and in extracellular medium) is an increasing function of IPTG induction. As our wetlab did not succeed in quantifying the secreted NGF, it is hard to figure out whether or not the desired concentration was obtained, but if our assumptions are valid, it could be reached with reasonable IPTG concentrations. Production of NGF with the tag has been detected by Mass spectrometry.

Our model is based on assumptions but it shows that within realistic parameters values, we can reasonably expect to obtain the optimal NGF concentration needed for neurite growth in the microfluidic chamber and it consequently paves the way to a functional proof of concept.

• M. Stamatakis and N. V. Mantzaris, “Comparison of deterministic and stochastic models of the lac operon genetic network,” Biophys. J., vol. 96, no. 3, pp. 887–906, 2009.
  
  
• A. Y. Weiße, D. A. Oyarzún, V. Danos, and P. S. Swain, “Mechanistic links between cellular trade-offs, gene expression, and growth,” Proc. Natl. Acad. Sci., vol. 112, no. 9, pp. E1038–E1047, 2015.
  
  
• R. Milo, “Useful fundamental BioNumbers handout.doc,” pp. 1–2, 2008.
  
  
• M. S. Packer, H. A. Rees, and D. R. Liu, “Phage-assisted continuous evolution of proteases with altered substrate specificity,” Nat. Commun., vol. 8, no. 1, 2017.
  
  
• H. Benabdelhak et al., “A specific interaction between the NBD of the ABC-transporter HlyB and a C-terminal fragment of its transport substrate haemolysin A,” J. Mol. Biol., vol. 327, no. 5, pp. 1169–1179, 2003.
  
  

Fick’s diffusion law

To simulate NGF diffusion in the microfluidic chip we consider a unidimensional conduit of axe x and a constant concentration of NGF introduced at one end of the canals. In this part, diffusion is assumed to be the only mechanism producing the gradient decay in the micro canals. We can model the diffusion characteristics of NGF with Fick’s second law of diffusion:
du / dt = C_diff d²u / dx²    (1)

Cdiff is assumed to be constant inside the conduit and depends on the material used.
There are also two boundary conditions:
at x=0:    du / dx |_(0,t)   (2)
at x=L:    du / dx |_(L,t)   (3)

Indeed, in the same material, the rate transfer of the diffusing NGF through the cross section of the micro canal is proportional to the concentration gradient normal to the cross section. It is assumed that the leakage of NGF at both ends of the micro canal is negligible because there should be little NGF at the ends the micro canals compared to the total amount of NGF and second because of a low NGF diffusion rate. The equation (1) can be solved with Euler’s method and we find the NGF concentration gradient at the position x and time t. The MatLab code is the following:

We obtain the following graphs:

Optimisation of the NGF gradient

Analysis of the model

To validate the model, we vary the three parameters (L, t_final, C_diff) to verify if the program corresponds to a diffusion phenomenon described in Fick’s second law of diffusion.

Observations:

1. When the length of the conduit increases but the duration of the experiment is fixed the NGF doesn’t have the time to diffuse in the entire conduit.
2. For instance, with a t_final= 3 600s the NGF molecules can’t diffuse further than x=0.2cm.

The higher the diffusion coefficient, the faster the molecules will diffuse in the conduit. Indeed, we observe in the model that with a fixed t_final:

1. NGF concentration at x=0.1 cm is 675 000 ng.ml^-1 for a diffusion coefficient C_diff = 15*10^-7 cm².s^-1
2. For a diffusion coefficient two times lower, the NGF concentration is 380 ng.ml¹

The results confirm the prediction of the Fick’s law model.

When the time length of the experiment lasts from 1 hour to 2 hours, the concentration of NGF is almost homogeneous in the entire conduit. At the end of the conduit, for x= 0.1 cm, the concentration of NGF equals to 910 ng.ml-1 when t_final= 7 200s whereas the concentration is 3 900 ng.ml^-1 when t_final=3 600s.

It is interesting to observe that when the duration of the experiment increases, the stationary regime is established: the NGF concentration in the conduit becomes independent of the position and time. Indeed, the concentation gradient of NGF in the conduit moves toward 0 for any position.

Explanation of the model

Baseline growth rate:

In our mathematical model, neurites grow at a constant growth rate defined as the baseline growth rate G0 when the concentration is below the threshold (assumed to be 995 ng.mL^-1). Neurites stop growing when the NGF concentration is higher than the threshold concentration. The value for the baseline growth rate G0 has been fixed at 20 μm.h^-1 for this model.

Concentration Gradient:

The extent of directional guidance is gradient steepness-dependent provided that the concentration gradient reaches the threshold value. The gradient factor k is a gradient steepness-dependent positive effect on the neurite growth rate.

In this model we assume that the baseline growth rate and the growth rate in the presence of concentration gradient follow an additive rule. This can be explained by the fact that both the NGF concentration and the its gradient can both individually contribute to neurite extension. The equation governing neurite outgrowth thus becomes:

dg / dt = G₀ + k ∂u / ∂x |_(g(t),t)       (4)

We can introduce a time parameter Tlag because the time taken to transmit the NGF signal is finite. The experiments show that the time lag for the cells to respond to NGF is approximately 1 day. The experiments show:
if t ≤ T_lag :     dg / dt = 0
else, if: t ≥ T_lag :     dg / dt = G₀ + k ∂u / ∂x |_(g(t),t)

Solving the model

To solve the equation (4) we are using Euler’s method forward because the gradient concentration of NGF depends on the length of the neurite (since neurites consume NGF).

The Equation (4):    
dg / dt = G₀ + k ∂u / ∂x |_(g(t),t)

Can be written as:    

g' = G₀ + k*f(g,t)

Which can be written as :    

g_n+1 - g_n / dt |_(g(t),t) = G_θ + k ∂u / ∂x |_(g(t),t)

We can therefore have an expression of g_n+1:    

g_n+1 = g_n + dt*[G₀ + k ∂u / ∂x |_(g(t),t)]

With initial values of g_θ, t_θ and ∂u / ∂x |_(g(t),t) we can find all the values of the g