Team:Pasteur Paris/Model

""

General introduction

The aim of our mathematical model is to simulate the growth of neurons towards our biofilm in response to the presence of pro Nerve Growth Factor (proNGF). proNGF is part of a family of proteins called neurotrophins. They are responsible for the development of new neurons, and for the growth and maintenance of mature ones. We created a deterministic model to help the wet lab establish the optimal concentration gradients of proNGF needed for the regrowth of the nerves. proNGF concentration and concentration gradient are key parameters affecting the growth rate and direction of neurites. Neurites growth has shown to be proNGF dose-dependent: if proNGF concentration is too low or too high, the growth rate is attenuated. In order to visualize the results of the model on a microfluidic chip, we used MATLAB, App Designer, Python, Gmsh, SpaceClaim and FreeFem. This is an important part of our project since it creates the link between the wet lab and dry lab.

We divided our model in three parts:

  1. Production of proNGF by the E. coli genetically modified
  2. Simulation of the diffusion of proNGF in a given environment
  3. Neurons growth in the presence of proNGF

Context of our model

Our project aims at creating a biofilm composed of genetically modified E. coli able to release a neurotrophic factor: proNGF. It helps to accelerate the connection between the neurons and the implant of the prosthesis; hence aiming at connecting the prosthesis and the amputee's neurons directly. This will enable the patient to have a more instinctive control of his prosthetic device. The nerves will be guided towards a conductive membrane surrounding our genetically modified biofilm. This membrane will then pass the neural signal of the regenerated nerves towards the electronic chip of the implant through wires. It will allow the patient to have a more instinctive and natural control than any other current prosthesis, and a reduced re-education time.

The aim of the wet lab is to test the biofilm on a microfluidic chip as a proof of concept. The chip is composed of two compartments: one contains the genetically modified E. coli that produce proNGF and the other one contains neurons. Microchannels link the two compartments in the middle of the chip, allowing the diffusion of proNGF and the growth of the neurites. Our model will hence be established on a microfluidic chip shape in order to share our results with the wet lab and indicate them the optimal concentration of proNGF needed according to our model.

We introduce different parameters in order to create our model :

g Length of the neurite outgrowth
dg/dt
 		Neurite outgrowth rate
u(x,t) Concentration of proNGF at the position x and time t
du/dt
 		proNGF concentration gradient at the position x and time t
Cdiff Diffusion coefficient of proNGF
K Gradient factor (growth rate of the neurite under the stimulation of the proNGF concentration gradient)
Gθ Baseline growth rate (neurite growth rate in absence of proNGF concentration gradient)
L Length of the conduit

proNGF production by genetically modified E. coli

As we want to obtain the best fitted proNGF concentration, we first simulate the production and secretion of our recombinant proNGF by transformed E. coli, in order to help the wetlab to optimize the induction and obtain the desired concentration, and to check whether we can theoretically obtain the optimal concentration for neurite growth.

Model Description

In this model, we include transcription, translation, translocation through E. coli membrane, protein folding and mRNA and protein degradation in cytoplasm and medium. proNGF synthesis is placed under Plac promoter, so we also modeled the IPTG induction. Finally, proNGF is secreted in the medium through Type I secretion system in which the export signal peptide is not cleaved during translocation. Our Biobrick is designed to synthetize and export TEV protease in order to cleave signal peptide and thus produce functional proNGF.

The molecular mechanism included in our model appears schematically in Figure 1.

Figure 1: Secretion mechanism of TEV and proNGF by our engineered bacteria

Our model includes the following variables:

Name Meaning
Iex IPTG outside the cell
Iin IPTG in the cytoplasm
Po Plac promoter occupied by repressor, prevent transcription
Pf Plac promoter with free lacO site
m mRNA for TEV and proNGF
m-r Ribosome-bound mRNA
proNGFc proNGF in cytoplasm
TEVc TEV protease in cytoplasm
(N-T)c proNGF-TEV complex in cytoplasm
proNGFcc Cleaved proNGF in cytoplasm, cannot be exported
proNGFt proNGF bound to transporter channel
TEVt TEV bound to transporter channel
t Transmembrane transporter
proNGFum Unfolded proNGF in medium with export peptide
proNGFm Folded proNGF in medium with export peptide
N-Tm Complex between proNGF with export peptide and functional TEV
TEVm TEV in medium with export peptide
proNGFf Functional proNGF in the medium

1. proNGF and TEV synthesis in the cytoplasm

The synthesis of proNGF and TEV is placed under the control of Plac promoter. The promoter can be in two different states: occupied (Po) by the repressor lacI, preventing RNA polymerase from binding and thus preventing transcription, or free (Pf) thanks to IPTG binding to the repressor. We assume that one IPTG molecule binds with one repressor molecule, 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 reactions with the same kinetic constant.

IPTG is not considered to be degraded neither in the cytoplasm nor in the medium.

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

For the TEV and proNGF translation, we first consider binding of ribosomes to ribosome binding site (the same association constant is used since the r.b.s. are the same), and then translation rate is proportional to the protein length. Since TEV and proNGF have approximately the same 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 proNGF has TEV cleaving site between the coding sequence and the export peptide, a fraction of proNGF 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, binds to its substrate and then cleaves the export peptide. This process can thus be modeled by the following equations:

K1, k-1 and k2 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.

2. proNGF and TEV secretion to the medium

The transport of proNGF 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, bounded 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 (proNGF and TEV) via its export signal peptide HlyA can bind to the HlyB-HlyD complex pore, forming a protein-transporter complex (proNGFt or TEVt). Translocation corresponds to the dissociation of this complex, resulting in restoring a free transporter and secreting proNGF or TEV in the medium (proNGFum and TEVm), which are the products.

3. Including growth rate

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 we need to put in our chip to produce the appropriate proNGF concentration.

4. proNGF folding and export peptide cleavage by TEV

Once in the medium, both proNGF and TEV are still bounded 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, bind to its substrate (which can be either proNGF with its export peptide or TEV with its export peptide) and then cleave the export peptide. This process can thus be modeled by the following equations:

5. mRNA and protein degradation

Finally, in cytoplasm and in the medium, mRNA and protein 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 (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 the following table.

NAME DESCRIPTION VALUE UNIT SOURCE
kt IPTG diffusion rate across the membrane 0.92 min-1 [1]
ki Association rate for derepression mechanism by IPTG 3 x 10-5 nM-1min-1 [1]
k-i Dissociation rate for derepression mechanism 4.8 x 103 min-1 [1]
α Transcription rate 2 mRNA.min-1nM-1 [3]
kr Association rate of ribosome with r.b.s 1 min-1mRNA-1 [2]
k-r Dissociation rate of ribosome with r.b.s 1 min-1 [2]
β Translation rate 4 nM.min-1mRNA-1 [3]
k1 Association rate of TEV with its substrate in the cytoplasm 7.8 x 10-7 min-1nM-1 Estimated from [4]
k-1 Dissociation rate of TEV with its substrate in the cytoplasm 6 x 10-4 min-1 Estimated from [4]
k2 Cleaving rate by TEV in cytoplasm 1.38 x 10-2 min-1 Estimated from [4]
k3 Association rate of proNGF and TEV with transmembrane transporter 6 x 10-4 min-1nM-1 [5]
k-3 Dissociation rate of proNGF and TEV with transporter 2.34 min-1 [5]
k4 Translocation rate within the transporter 2.1 min-1 [5]
kf proNGF folding rate in the medium 0.28 min-1
k5 Association rate of TEV with its substrate in the medium 7.8 x 10-5 min-1nM-1 [4]
k-5 Dissociation rate of TEV with its substrate in the medium 0.06 min-1nM [4]
k6 Cleaving rate by TEV in the medium 1.38 min-1nM-1 [4]
δm mRNA degradation rate 0.462 min-1 [1]
δpc Protein degradation rate in cytoplasm 0.2 min-1 [1]
δpm Protein degradation rate in extracelular medium 0.1 min-1 [1]

MODEL RESULTS

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

Figure 2: Comparison of cytoplasmic and secreted NGF with a single-cell model (IPTG induction 1 mM)

After the initial dynamics, concentration of secreted proNGF quickly reaches a steady state , which is then only driven by the bacterial population dynamics. If we consider a bacterial culture in stationary phase, we can consequently consider that the initial proNGF concentration is constant. Our model predicts that the majority of recombinant protein remains cytoplasmic or is secreted but not functional (we consider as "non-functional proNGF" the recombinant proteins that are not folded or still have a C-terminal HlyA signal peptide), as it appears in Fig1.

The aim of this first model is to demonstrate that we can expect an appropriate secreted recombinant proNGF 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 proNGF 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 proNGF secretion with what we need.

Influence of number of transporters

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 proNGF concentration we can expect.

The following graph shows the predicted proNGF 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.):

Figure 3: Comparison of cytoplasmic and secreted NGF when the number of transporters varies

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 proNGF concentration we can expect.

The following graph shows the predicted proNGF 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.):

Influence of translocation rate

Figure 4: Secreted NGF as a function of translocation rate

As expected, the more transporters the cell has, the more recombinant proNGF is secreted, but the amount of functional secreted proNGF (in blue) remains limited due to TEV protease cleaving efficiency.

Taking in account the number of E. coli cells and the dilution factor between intracellular and extracellular space, we obtain for 500 transporters a concentration of functional proNGF of 1 nM, which corresponds to 24 ng/mL. This is still 10 times lower than what we need to observe neurite growth. Enhancing signal peptide cleavage by a more efficient enzyme should help solve the problem, since we could expect 5 nM functional proNGF if the totality of the secreted proNGF was cleaved.

IPTG induction level

One of the parameters our wetlab team is able to adjust is IPTG induction in the microchannel chip in order to optimize the obtained proNGF concentration. Consequently, we studied the dependence of secreted proNGF with IPTG initial concentration.

As expected the final proNGF 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 proNGF, 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 proNGF with the tag has been detected by Mass spectrometry.

Figure 5: Comparison of cytoplasmic and secreted NGF for different IPTG induction level

PERSPECTIVES

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

Next modeling steps:

  • It would be worth isolating and quantifying secreted recombinant proNGF in order to confront model and experiments, and be able to determine some of the kinetics parameters values we used (such as translocation rate)
  • This program is designed to model the microchip proof-of-concept experiment but we will adapt it to our final biofilm device to predict its behavior

proNGF diffusion diffusion in a given environment


We are looking to understand the way the proNGF spreads inside the conduit once it is produced. This will help us to determine the proNGF concentration u(x,t) (ng.mL-1) as a function of the distance x (cm) from the production site of proNGF.

Fick’s diffusion law

To simulate proNGF diffusion in the microfluidic chip we consider a unidimensional conduit of axe x and a constant concentration of proNGF 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 proNGF with Fick’s second law of diffusion:
du / dt = Cdiff d2u / dx2    (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 proNGF 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 proNGF at both ends of the micro canal is negligible because there should be little proNGF at the ends the micro canals compared to the total amount of proNGF and second because of a low proNGF diffusion rate. The equation (1) can be solved with Euler’s method and we find the proNGF concentration gradient at the position x and time t. The MatLab code is the following:

We displayed our results showing a decrease of the concentration of proNGF (u(x,t)) depending on the distance of the conduit x.

We used the following parameters for the model:

Length of the conduit: L 0.1 cm
Diffusion coefficient of proNGF : Cdiff 7,8*10-7 cm2.s-1
Time length of the experiment: t_final 3 600 s

We obtain the following graphs:

Optimisation of the proNGF gradient

To optimize the accuracy of the proNGF gradient we interpolate the curve u(x)=f(x). Consequently, we obtain the f polynomial function easier to derive and a polynomial function of the gradient with a better accuracy than with the first method. The program is the following:

With the same parameters as with the previous model we obtain the following graphs:

Analysis of the model

To validate the model, we vary the three parameters (L, t_final, Cdiff) 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 proNGF doesn’t have the time to diffuse in the entire conduit.
  2. For instance, with a t_final= 3 600s the proNGF 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. proNGF concentration at x=0.1 cm is 675 000 ng.ml-1 for a diffusion coefficient Cdiff = 15*10-7 cm2.s-1
  2. For a diffusion coefficient two times lower, the proNGF concentration is 380 ng.ml1

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 proNGF is almost homogeneous in the entire conduit. At the end of the conduit, for x= 0.1 cm, the concentration of proNGF 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 proNGF concentration in the conduit becomes independent of the position and time. Indeed, the concentation gradient of proNGF in the conduit moves toward 0 for any position.

Neurons growth in the presence of proNGF


In this part our goal is to determine the length of the neurite outgrowth (g(t)) in response to the gradient concentration of proNGF.

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 proNGF 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 proNGF concentration and the its gradient can both individually contribute to neurite extension. The equation governing neurite outgrowth thus becomes:

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

We can introduce a time parameter Tlag because the time taken to transmit the proNGF signal is finite. The experiments show that the time lag for the cells to respond to proNGF is approximately 1 day. The experiments show:
if t ≤ Tlag :     dg / dt = 0
else, if: t ≥ Tlag :     dg / dt = G0 + 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 proNGF depends on the length of the neurite (since neurites consume proNGF).

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

Can be written as:    

g' = G0 + k*f(g,t)

Which can be written as :    

gn+1 - gn / dt |(g(t),t) = Gθ + k ∂u / ∂x |(g(t),t)

We can therefore have an expression of gn+1:    

gn+1 = gn + dt*[G0 + 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

REFERENCES

  • 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.

  • Defining the concentration gradient of nerve growth factor for guided neurite outgrowth, XCao M.SShoichet, March 2001

  • Immobilized Concentration Gradients of Neurotrophic Factors Guide Neurite Outgrowth of Primary Neurons in Macroporous Scaffolds, Moore K, MacSween M, Shoichet M, feb 2006

  • Mathematical Modeling of Guided Neurite Extension in an Engineered Conduit with Multiple Concentration Gradients of Nerve Growth Factor (proNGF), Tse TH, Chan BP, Chan CM, Lam J, sep 2007

  • Mathematical modeling of multispecies biofilms for wastewater treatment, Maria Rosaria Mattei, november 2005