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.

Our model includes the following variables:

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.

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.

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.

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:

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.

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.

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.

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.):

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.):

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.

Fick’s diffusion law

To simulate proNGF diffusion in the microfluidic chip we consider a unidimensional conduit of axe x (cm) and a constant concentration rate 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. This corresponds to a unidimensional version of equation (1) with s = 0:
du / dt = C_diff d²u / dx²    (1)

C_diff is assumed to be constant inside the conduit and depends on the material used.

The Neumann boundary conditions of this equation are:
at x=0:    du / dx |_(0,t) = u₀   (2)
at x=L:    du / dx |_(L,t) = 0   (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 obtain the following graphs:

Optimisation of the proNGF gradient

Analysis of the model

To validate the model, we vary two parameters (L and 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 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 C_diff = 15*10^-7 cm².s^-1
2. For a diffusion coefficient two times lower, the proNGF 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 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.

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 its gradient can 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 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 ≤ T_lag :     dg / dt = 0
else, if: t ≥ T_lag :     dg / dt = G₀ + k ∂u / ∂x |_(g(t),t)

Euler's Method

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

Solving the Model

We noticed that there was a mistake in the article regarding the gradient steepness-dependent factor k. Therefore, by examining the results of the article, we were able to find a new coherent value of k, different from what was originally written in the article. To see the details of the mistake and our strategy to find the right value, click here.

In the following graphs, the red curve corresponds to a nerve growth inside a unidimensionnal canal without any NGF while the blue one corresponds to the situation where there is an NGF gradient inside the canal.

Our strategy to compute this phenomenon is the following.

For each different time, (spaced by the value dt) the position (in cm) of an axon is put inside two unidimensional matrices, g and g_control. The matrix g holds the values of positions when there is a gradient of NGF, while there isn’t for g_control. At each time and for the corresponding position, the script we used to calculate the gradient of NGF (in part 2 ) is ran with those new parameters. As we obtain the value of the gradient of NGF at this time and at a position g_n, we can calculate the new position at the end of the axon g_n+1 by using the formula written above.

The first set of parameters we use is the following :

Not only does the result clearly demonstrates the need of an NGF gradient in order to boost the nerve growth, it also indicates how fast will a neuron grow in a certain amount of time, knowing the NGF concentration at the beginning of the canal.

We then carried on the study regarding two different parameters.

We studied the influence of the length of the canal on the nerve growth. We worked with 2 sets of parameters:

The more the length of the device rises, the longer it will take for the NGF concentration to be homogenous inside the canal. Therefore, to boost the growth nerve by having a gradient concentration of NGF in the media, the length of the device must be taken into account. The difference is significant since, at t=50 000 s, the magnitude of the gradient concentration of NGF in the canal worth 10³ when L = 0.4 cm while it worth 10^-5 when L= 0.1 cm.

If the length is too small, the gradient concentration won’t be high enough to boost the growth nerves significantly, as shown in the graph Evolution of the nerve growth above.

We studied the influence of the length of the canal on the nerve growth. We worked with two sets of parameters

The results are shown below.

The results show that, the more the initial concentration of NGF increases, the more its gradient reaches higher values and therefore the faster the neurons grow. It would appear that increasing the initial concentration of NGF would help to boost the nerve growth. However, if the concentration of NGF is too high, it would cause the opposite effect as neurons would start to die. Finding the right compromise, depending on the length of our final device and the coefficient diffusion of NGF in the future media, will be of paramount importance. Finding this compromise will be possible thanks to our model.

Our modeling work has permitted us to study the secretion, diffusion and influence of NGF on the growth of the neurons. With our model, we are able to optimize two parameters : the length of the micro channels and the initial concentration. The wet lab took into consideration our results to do the experiments on the influence of the concentration of NGF on the growth of the axons. Indeed, they observed the growth of the axons of E18 cortex cells for different concentrations of NGF: 0, 50, 250, 500, 750 and 900 ng/mL. The wet lab’s results were coherent with our model.

The experiments show that until a certain concentration the growth of the neurons increases with the presence of NGF. For a concentration between 250 and 750 ng/mL the presence of NGF increases significantly the growth of the axons. On the model (figure 9 and 10), the higher the NGF concentration, the higher the gradient concentration of NGF so the faster the neurons will grow.

Figure 9: Experimental results of the wet lab showing the importance of the growth of the axons for different concentrations of NGF. β-III Tubulin is a bio marker of neuron cell differentiation that indicates the growth of axons

Figure 10: The modeling results showing the nerves growth (cm) for different concentrations of NGF (purple: 0 ng/mL, blue: 250 ng/mL, red: 500 ng/mL, yellow: 750 ng/mL)

The model helped the wet lab establish the concentration limit of NGF above which the NGF doesn’t have any more influence on the growth of the neurons. The wet lab’s concentration limit is coherent with ours: their concentration limit is approximatevely 900 ng/mL whilst the model shows a concentration limit of 995 ng/mL.

The wet lab has done the series of experiments on a 96 wells plate in order to optimize the number of samples. The next step for the wetlab is to experimentally verify the influence of the length of the microchannels in the microfluidic chip on the growth of the nerves. The model is able to provide information on the optimization of the length of the microchannels which could be of use for the wet lab. Another improvement would be to calculate the diffusion coefficient in the microfluidic chip media.

Our model will be used to prototype the final device to help and establish the NGF concentration needed to control nerves’ growth. The length of the nerves needed to reach the interface depends on individual. As mentioned in the design scenario, chemical induction for bacteria regarding NGF production might be considered. Since the model manages to link induction to diffusion to nerves growth, it will enable to know how much NGF needs to be produced for each individual.

The next step consists in keep trying to get in touch with the authors of the article or contacting other experts to make our model completely fulfill its major role in NeuronArch.