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