Team:Valencia Biocampus/Modeling

Biobricks modeling

In this section, clumping and PHA production are modeled and fitted to a equation. These things will be considered in other sections.

Clumping

PHA production

C. elegans behavior

In order to understand the system as a whole, we started with the simplest case, so our first study is based on the behavior of a single worm, initially in the absence of chemoattractant and later adding it.

Our first approach consists in the study of a single worm in a surface without chemoattractant, under this condition it is observed that a C. elegans moves randomly, this is defined as a Random Walk.

Random Walk is the mathematical formalization of a trajectory that consists of taking successive random steps. We ran several simulations in Scilab and C++, in order to represent a single C. elegans moving this way in the absence of attractants.

Random walks are used to describe the trajectories of many motile animals and microorganisms. They are useful for both qualitative and quantitative descriptions of the behavior of such creatures. In our case a single C. elegans was simulated as a point: $(x_t, y_t)$
At each instant $t$ in the simulation, step length $l_t$ and direction $\theta_t$ the equations are the following:
$$l_t = \Delta t\;\nu_t$$ $$\theta_t = \theta_{t-1} + \Delta t \left(\frac{d\theta_t}{dt} + \delta\right) $$ Where $\Delta t$ is the duration of the time step, $\nu_t$ is the instantaneous speed, $\frac{d\theta_t}{dt}$ is the instantaneous turning rate (using the convention that $\frac{d\theta_t}{dt}$$ > 0$ is a right turn and $\frac{d\theta_t}{dt}$$ < 0$ is a left turn) and $\delta$ is the turning bias. The values of these variables were taken randomly from different gaussian distributions, that we identified by sampling but also obtained from some papers:

$\nu_t$: normal distribution with $\sigma$ = 0.0152 cm/s and $\mu$ = 0.00702 cm/s
\dot{$\theta$}: normal distribution with $\sigma$ = 0.0150273 rad/s and $\mu$ = 0.6789331 rad/s
$\tdelta$: normal distribution with $\sigma$ = 0.0076969 rad/s and $\mu$ = 0.0370010 rad/s

Simulations were performed in Scilab and C++, proceeding as described in (Pierce-Shimomura et al., 1999). We contacted the authors but they could not provide us with the code, so we developed our own script based on the article, obtaining similar results:

Add simulation pictures

Chemotaxis of a single worm

However, we are interested in the movement of the nematode in a gradient of chemoattractant. Studies revealed that the behavior of E. coli during chemotaxis is remarkably similar to the behavior of a single C. elegans in the presence of a chemotactic source (Bargmann, 2006), a mechanism called the pirouette model in C. elegans (Pierce-Shimomura et al., 1999) and the biased random walk in bacteria (Berg, 1993; Berg, 1975). The basis of these models is a strategy that uses a short-term memory of attractant concentration to decide whether to maintain the current direction of movement or to change to a new one.
Following this model, we defined the turning rate as a function of $\frac{\partial [C]}{\partial t}$. Where $[C]$ is the concentration ($[\;]$) of chemoattracant ($C$). Simulations were performed as expected, showing a bias in the random walk:

Add simulation pictures

Chemotaxis of a single worm consuming attractant

As a last step of the study of a single worm, we decided to modelize the most realistic behavior. This behavior happens when, in addition to moving randomly but biased toward the gradient of chemoattractant, we consider the fact that it is constantly consuming its “food”, so these gradients will not be any more constants in time.

In order to implement this model, we realized that was impossible to study still considering the space as an attractant function, cause it changes its shape in every time step, and it made simulations really slow (maybe, days of computation). So, we meshed the space, and gave every point a different weight, depending on time: that was the starting of the matrix and approximated (numerical) calculations. What we obtained was close to have relevant impacts in the course of our project, because we could interfere in the C. elegans path, and not only in its final position, making it moving from one source of chemoattractant to another:

Add simulation pictures
aquí va el vídeo
Moreover, these numerical methods opened ours eyes to study their movement as a whole, with partial differential equations.

Bargmann CI (2006) Chemosensation in C. elegans. Wormbook
Berg HC (1993) Random walks in biology. Princeton, NJ: Princeton UP.
Pierce-Shimomura JT, Morse TM, Lockery SR (1999) The Fundamental Role of Pirouettes in C. elegans Chemotaxis. The journal of Neuroscience, 19(21):9557-9569

Group behavior

In practice, a single C. elegans is not employed to perform the task but a group of worms. In this case, the random walk equations obtained for a single worm behavior can be transformed into partial differential equations (PDEs) to depict the distribution of the population, as the ones describing a diffusion process.
Mathematical demonstration of macroscopic RW as a diffusion process. Link?
The equations are the following: $$\frac{\partial [C]}{\partial t} = D\;\nabla^2[C] - \nabla\cdot(\overrightarrow{\nu}\;[C]) $$ For the 2-dimensional case, $[C] = [C(x, y, t)]$ is the concentration of C. elegans at a given instant $t$ and at a given point $(x, y)$, $D$ is the diffusion coefficient and $\overrightarrow{\nu}$ is the attraction field for the C. elegans, which basically stands for its velocity. This expression naturally expands to: $$\frac{\partial [C]}{\partial t} = D \; \left(\frac{\partial^2 [C]}{\partial x^2} + \frac{\partial^2 [C]}{\partial y^2}\right) - \left(\frac{\partial\nu_x}{dx}[C] + \nu_x\frac{\partial [C]}{\partial x} + \frac{\partial\nu_y}{dy}[C] + \nu_y\frac{\partial [C]}{\partial y}\right) $$
However, in our case we face a diffusion system with drift, in which the last term arises from the biased movement of C. elegans. To obtain a PDE that reflects this behavior, we employed difference equations using Taylor series and apropiate limits
Mathematical demonstration of C. elegans chemotactic behavior using PDEs from biased random walk employing the pirouette model. Link? Properties as mean or SD calculated from the equation?

PDEs. Parameters, as diffusion coefficient (D) or drift force (b) can be calculated from (Pierce-Shimomura, 1999)

We first started by building up an explicit finite difference method, this is numerically fast but it is prone to instabilities in the solutions. Therefore, we decided to develop a Crank-Nicolson method which is unconditionally stable, thus being suitable to carry out parameter identification of a model. Nevertheless, it has the drawback of being numerically more intensive, because a set of algebraic equations must be solved in each iteration.

The Crank-Nicolson method is obtained by using a forward difference in time ($t_n$) and central differences in space $(x_j, y_k)$. However, by giving discrete values for time and space $(t_n, x_j, y_k)$ to our Crank-Nicolson method and zeroing boundary conditions (since we have to define some limit for the plate) and also setting initial conditions to a specific initial distribution, we achieved the matrix recursive formula: $$[\mathcal{C}^{n+1}_{j, k}] = \mathcal{M}\;[\mathcal{C}^n_{j, k}]$$ Where $\mathcal{M}$ is a block tridiagonal matrix depending on the parameters and the velocities in $x$ and $y$ directions at some points in space.

Now, we can predict the behavior of a set of worms in presence of a chemoatractant. Simulations were run using Scilab. As seen, the distribution...

Add simulation

Wormboys in action

Once defined the behavior of a set of worms with no bacteria attached, we can model the action of microorganisms while they are riding C. elegans. In our model, bacterial substrate and worm chemoatractant are the same molecule, so it is consume while C. elegans is moving through the place. This alters the pattern of the molecule of interest, modifying the trajectory of the nematode. To take into account this feature of the system, a reaction term is added in the previous equation.

PDE

System performance

The main objective of the project is employing the equations to describe a new, theoretical, solid-batch biorreactor and comparing it versus liquid-batch biorreactor.