Team:Valencia Biocampus/Modeling
From 2013.igem.org
Modeling
The main goal of our modeling project is to accurately predict the behavior of our system in several issues, from the mounting of bacteria on C. elegans to the performance of our worms reaching the place of interest. In order to do that, we use several modeling techniques.
The movement of C. elegans in the presence of a chemoattractant in order to carry our bacteria to that source is the main issue of our project, so we consider modeling this aspect and employing it as scaffold for the whole modeling project. Several layers show up that it must be modeled in different ways. In our approach, we mathematically describe each layer, from the simplest to the most complex, integrating each one.
The workflow in the C. elegans movement is the following:
Single Worm Behavior | Single Worm Chemotaxis | Group Behavior | Wormboys in action | System performance | ||
---|---|---|---|---|---|---|
Considerations | Chemoattractant | ![]() |
![]() |
![]() |
![]() |
![]() |
Bacteria | ![]() |
![]() |
![]() |
![]() |
![]() |
|
Modeling approach | Random Walk | Biased Random Walk | Partial Differential Equations | Partial Differential Equations | Ordinary Differential Equations |
The model is improved adding the data obtained from the experiments in order to achieve a holistic model which can predict the distribution of our worms, the kinetics of the present bacteria and the concentration of the substrate in each moment.
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} (= \dot{\theta_t})$ 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_t}$: normal distribution with $\sigma = 0.0150273\;rad/s$ and $\mu = 0.6789331\;rad/s$
$\delta$: 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:
Simulated with our own software SimuElegans.
Please note that the actual movements of C. elegans are much slower and thus the simulation in this video is accelerated for convenience.
In these simulations we considered $\Delta t = 1\;s$, but actually simulated it every $0.01\;s$, so its $\frac{1\;s}{0.01\;s} = 100$ times accelerated.
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
Simulated with our own software SimuElegans.
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
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.
<a href= "https://2013.igem.org/Team:Valencia_Biocampus/Demonstration/Diffusion1">Mathematical demonstration of macroscopic RW as a diffusion process</a>
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_k$) and central differences in space $(x_i, y_j)$.
However, by giving discrete values for time and space $(t_k, x_i, y_j)$ 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}^{k+1}] = \mathcal{M}\;[\mathcal{C}^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
Description of the numerical method
1) The first step is to discretize our differential equation, replacing the time derivatives by the approximate forward differences and the space derivatives by central differences, being $\;q\;$ and $\;h\;$ the time and space step, respectively:
$$\; 2 \; \frac{C^{k+1}_{i,j} - C^{k}_{i,j}}{q}\;= \\ \;D \; \left(\frac{C^{k+1}_{i+1,j} - 2C^{k+1}_{i,j}+C^{k+1}_{i-1,j}}{h^2} + \frac{C^{k+1}_{i,j+1} - 2C^{k+1}_{i,j}+C^{k+1}_{i,j-1}}{h^2}\right) - \left(\frac{\nu_{x_{i+1,j}} - \nu_{x_{i-1,j}}}{2h} \; C^{k+1}_{i,j} + \nu_{x_{i,j}} \; \frac{C^{k+1}_{i+1,j} - C^{k+1}_{i-1,j}}{2h} + \frac{\nu_{y_{i,j+1}} - \nu_{y_{i,j-1}}}{2h} \; C^{k+1}_{i,j} + \nu_{y_{i,j}} \; \frac{C^{k+1}_{i,j+1} - C^{k+1}_{i,j-1}}{2h}\right) \\ + D \; \left(\frac{C^{k}_{i+1,j} - 2C^{k}_{i,j}+C^{k}_{i-1,j}}{h^2} + \frac{C^{k}_{i,j+1} - 2C^{k}_{i,j}+C^{k}_{i,j-1}}{h^2}\right) - \left(\frac{\nu_{x_{i+1,j}} - \nu_{x_{i-1,j}}}{2h} \; C^{k}_{i,j} + \nu_{x_{i,j}} \; \frac{C^{k}_{i+1,j} - C^{k}_{i-1,j}}{2h} + \frac{\nu_{y_{i,j+1}} - \nu_{y_{i,j-1}}}{2h} \; C^{k}_{i,j} + \nu_{y_{i,j}} \; \frac{C^{k}_{i,j+1} - C^{k}_{i,j-1}}{2h}\right)$$
Then, placing all (k+1) terms at left and (k) terms at right and recalling the constants $\;\lambda \; = \; \frac{q \; D}{2 \; h^2}\;$ and $\;\mu \; = \; \frac{q}{4 \; h}\;$ it results a more handy equation:
$$ C^{k+1}_{i+1,j} \; \left(- \lambda + \nu_{x_{i,j}} \; \mu \right) + C^{k+1}_{i,j+1} \; \left(- \lambda + \nu_{y_{i,j}} \; \mu \right) + C^{k+1}_{i,j} \; \left(1 + 4 \; \lambda + \nu^{(1,0)}_{{x}_{i,j}} + \nu^{(0,1)}_{{y}_{i,j}}\right) + C^{k+1}_{i-1,j} \; \left(- \lambda - \nu_{x_{i,j}} \; \mu \right) + C^{k+1}_{i,j-1} \; \left(- \lambda - \nu_{y_{i,j}} \; \mu \right) \; = \\ \; C^{k}_{i+1,j} \; \left(\lambda - \nu_{x_{i,j}} \; \mu \right) + C^{k}_{i,j+1} \; \left(\lambda - \nu_{y_{i,j}} \; \mu \right) + C^{k}_{i,j} \; \left(1 - 4 \; \lambda - \nu^{(1,0)}_{{x}_{i,j}} - \nu^{(0,1)}_{{y}_{i,j}}\right) + C^{k}_{i-1,j} \; \left(\lambda + \nu_{x_{i,j}} \; \mu \right) + C^{k}_{i,j-1} \; \left(\lambda + \nu_{y_{i,j}} \; \mu \right) $$
2) For being able to compute the method in a matrix way, we were forced to do a small cencrete example, in order to extrapolate then the procedure into a general method. For that, we meshed the space in three steps each coordinate, so that, i=1,2,3 and j=1,2,3:
$ i=1,\;j=1 $
$$ C^{k+1}_{2,1} \; \left(- \lambda + \nu_{x_{1,1}} \; \mu \right) + C^{k+1}_{1,2} \; \left(- \lambda + \nu_{y_{1,1}} \; \mu \right) + C^{k+1}_{1,1} \; \left(1 + 4 \; \lambda + \nu^{(1,0)}_{{x}_{1,1}} + \nu^{(0,1)}_{{y}_{1,1}}\right) + C^{k+1}_{0,1} \; \left(- \lambda - \nu_{x_{1,1}} \; \mu \right) + C^{k+1}_{1,0} \; \left(- \lambda - \nu_{y_{1,1}} \; \mu \right) \; = \\ \; C^{k}_{2,1} \; \left(\lambda - \nu_{x_{1,1}} \; \mu \right) + C^{k}_{1,2} \; \left(\lambda - \nu_{y_{1,1}} \; \mu \right) + C^{k}_{1,1} \; \left(1 - 4 \; \lambda - \nu^{(1,0)}_{{x}_{1,1}} - \nu^{(0,1)}_{{y}_{1,1}}\right) + C^{k}_{0,1} \; \left(\lambda + \nu_{x_{1,1}} \; \mu \right) + C^{k}_{1,0} \; \left(\lambda + \nu_{y_{1,1}} \; \mu \right) $$
$ i=1,\;j=2 $
$$ C^{k+1}_{2,2} \; \left(- \lambda + \nu_{x_{1,2}} \; \mu \right) + C^{k+1}_{1,3} \; \left(- \lambda + \nu_{y_{1,2}} \; \mu \right) + C^{k+1}_{1,2} \; \left(1 + 4 \; \lambda + \nu^{(1,0)}_{{x}_{1,2}} + \nu^{(0,1)}_{{y}_{1,2}}\right) + C^{k+1}_{0,2} \; \left(- \lambda - \nu_{x_{1,2}} \; \mu \right) + C^{k+1}_{1,1} \; \left(- \lambda - \nu_{y_{1,2}} \; \mu \right) \; = \\ \; C^{k}_{2,2} \; \left(\lambda - \nu_{x_{1,2}} \; \mu \right) + C^{k}_{1,3} \; \left(\lambda - \nu_{y_{1,2}} \; \mu \right) + C^{k}_{1,2} \; \left(1 - 4 \; \lambda - \nu^{(1,0)}_{{x}_{1,2}} - \nu^{(0,1)}_{{y}_{1,2}}\right) + C^{k}_{0,2} \; \left(\lambda + \nu_{x_{1,2}} \; \mu \right) + C^{k}_{1,1} \; \left(\lambda + \nu_{y_{1,2}} \; \mu \right) $$
$ i=1,\;j=3 $
$$ C^{k+1}_{2,3} \; \left(- \lambda + \nu_{x_{1,3}} \; \mu \right) + C^{k+1}_{1,4} \; \left(- \lambda + \nu_{y_{1,3}} \; \mu \right) + C^{k+1}_{1,3} \; \left(1 + 4 \; \lambda + \nu^{(1,0)}_{{x}_{1,3}} + \nu^{(0,1)}_{{y}_{1,3}}\right) + C^{k+1}_{0,3} \; \left(- \lambda - \nu_{x_{1,3}} \; \mu \right) + C^{k+1}_{1,2} \; \left(- \lambda - \nu_{y_{1,3}} \; \mu \right) \; = \\ \; C^{k}_{2,3} \; \left(\lambda - \nu_{x_{1,3}} \; \mu \right) + C^{k}_{1,4} \; \left(\lambda - \nu_{y_{1,3}} \; \mu \right) + C^{k}_{1,3} \; \left(1 - 4 \; \lambda - \nu^{(1,0)}_{{x}_{1,3}} - \nu^{(0,1)}_{{y}_{1,3}}\right) + C^{k}_{0,3} \; \left(\lambda + \nu_{x_{1,3}} \; \mu \right) + C^{k}_{1,2} \; \left(\lambda + \nu_{y_{1,3}} \; \mu \right) $$
$ i=2,\;j=1 $
$$ C^{k+1}_{3,1} \; \left(- \lambda + \nu_{x_{2,1}} \; \mu \right) + C^{k+1}_{2,2} \; \left(- \lambda + \nu_{y_{2,1}} \; \mu \right) + C^{k+1}_{2,1} \; \left(1 + 4 \; \lambda + \nu^{(1,0)}_{{x}_{2,1}} + \nu^{(0,1)}_{{y}_{2,1}}\right) + C^{k+1}_{1,1} \; \left(- \lambda - \nu_{x_{2,1}} \; \mu \right) + C^{k+1}_{2,0} \; \left(- \lambda - \nu_{y_{2,1}} \; \mu \right) \; = \\ \; C^{k}_{3,1} \; \left(\lambda - \nu_{x_{2,1}} \; \mu \right) + C^{k}_{2,2} \; \left(\lambda - \nu_{y_{2,1}} \; \mu \right) + C^{k}_{2,1} \; \left(1 - 4 \; \lambda - \nu^{(1,0)}_{{x}_{2,1}} - \nu^{(0,1)}_{{y}_{2,1}}\right) + C^{k}_{1,1} \; \left(\lambda + \nu_{x_{2,1}} \; \mu \right) + C^{k}_{2,0} \; \left(\lambda + \nu_{y_{2,1}} \; \mu \right) $$
$ i=2,\;j=2 $
$$ C^{k+1}_{3,2} \; \left(- \lambda + \nu_{x_{2,2}} \; \mu \right) + C^{k+1}_{2,3} \; \left(- \lambda + \nu_{y_{2,2}} \; \mu \right) + C^{k+1}_{2,2} \; \left(1 + 4 \; \lambda + \nu^{(1,0)}_{{x}_{2,2}} + \nu^{(0,1)}_{{y}_{2,2}}\right) + C^{k+1}_{1,2} \; \left(- \lambda - \nu_{x_{2,2}} \; \mu \right) + C^{k+1}_{2,1} \; \left(- \lambda - \nu_{y_{2,2}} \; \mu \right) \; = \\ \; C^{k}_{3,2} \; \left(\lambda - \nu_{x_{2,2}} \; \mu \right) + C^{k}_{2,3} \; \left(\lambda - \nu_{y_{2,2}} \; \mu \right) + C^{k}_{2,2} \; \left(1 - 4 \; \lambda - \nu^{(1,0)}_{{x}_{2,2}} - \nu^{(0,1)}_{{y}_{2,2}}\right) + C^{k}_{1,2} \; \left(\lambda + \nu_{x_{2,2}} \; \mu \right) + C^{k}_{2,1} \; \left(\lambda + \nu_{y_{2,2}} \; \mu \right) $$
$ i=2,\;j=3 $
$$ C^{k+1}_{3,3} \; \left(- \lambda + \nu_{x_{2,3}} \; \mu \right) + C^{k+1}_{2,4} \; \left(- \lambda + \nu_{y_{2,3}} \; \mu \right) + C^{k+1}_{2,3} \; \left(1 + 4 \; \lambda + \nu^{(1,0)}_{{x}_{2,3}} + \nu^{(0,1)}_{{y}_{2,3}}\right) + C^{k+1}_{1,3} \; \left(- \lambda - \nu_{x_{2,3}} \; \mu \right) + C^{k+1}_{2,2} \; \left(- \lambda - \nu_{y_{2,3}} \; \mu \right) \; = \\ \; C^{k}_{3,3} \; \left(\lambda - \nu_{x_{2,3}} \; \mu \right) + C^{k}_{2,4} \; \left(\lambda - \nu_{y_{2,3}} \; \mu \right) + C^{k}_{2,3} \; \left(1 - 4 \; \lambda - \nu^{(1,0)}_{{x}_{2,3}} - \nu^{(0,1)}_{{y}_{2,3}}\right) + C^{k}_{1,3} \; \left(\lambda + \nu_{x_{2,3}} \; \mu \right) + C^{k}_{2,2} \; \left(\lambda + \nu_{y_{2,3}} \; \mu \right) $$
$ i=3,\;j=1 $
$$ C^{k+1}_{4,1} \; \left(- \lambda + \nu_{x_{3,1}} \; \mu \right) + C^{k+1}_{3,2} \; \left(- \lambda + \nu_{y_{3,1}} \; \mu \right) + C^{k+1}_{3,1} \; \left(1 + 4 \; \lambda + \nu^{(1,0)}_{{x}_{3,1}} + \nu^{(0,1)}_{{y}_{3,1}}\right) + C^{k+1}_{2,1} \; \left(- \lambda - \nu_{x_{3,1}} \; \mu \right) + C^{k+1}_{3,0} \; \left(- \lambda - \nu_{y_{3,1}} \; \mu \right) \; = \\ \; C^{k}_{4,1} \; \left(\lambda - \nu_{x_{3,1}} \; \mu \right) + C^{k}_{3,2} \; \left(\lambda - \nu_{y_{3,1}} \; \mu \right) + C^{k}_{3,1} \; \left(1 - 4 \; \lambda - \nu^{(1,0)}_{{x}_{3,1}} - \nu^{(0,1)}_{{y}_{3,1}}\right) + C^{k}_{2,1} \; \left(\lambda + \nu_{x_{3,1}} \; \mu \right) + C^{k}_{3,0} \; \left(\lambda + \nu_{y_{3,1}} \; \mu \right) $$
$ i=3,\;j=2 $
$$ C^{k+1}_{4,2} \; \left(- \lambda + \nu_{x_{3,2}} \; \mu \right) + C^{k+1}_{3,3} \; \left(- \lambda + \nu_{y_{3,2}} \; \mu \right) + C^{k+1}_{3,2} \; \left(1 + 4 \; \lambda + \nu^{(1,0)}_{{x}_{3,2}} + \nu^{(0,1)}_{{y}_{3,2}}\right) + C^{k+1}_{2,2} \; \left(- \lambda - \nu_{x_{3,2}} \; \mu \right) + C^{k+1}_{3,1} \; \left(- \lambda - \nu_{y_{3,2}} \; \mu \right) \; = \\ \; C^{k}_{4,2} \; \left(\lambda - \nu_{x_{3,2}} \; \mu \right) + C^{k}_{3,3} \; \left(\lambda - \nu_{y_{3,2}} \; \mu \right) + C^{k}_{3,2} \; \left(1 - 4 \; \lambda - \nu^{(1,0)}_{{x}_{3,2}} - \nu^{(0,1)}_{{y}_{3,2}}\right) + C^{k}_{2,2} \; \left(\lambda + \nu_{x_{3,2}} \; \mu \right) + C^{k}_{3,1} \; \left(\lambda + \nu_{y_{3,2}} \; \mu \right) $$
$ i=3,\;j=3 $
$$ C^{k+1}_{4,3} \; \left(- \lambda + \nu_{x_{3,3}} \; \mu \right) + C^{k+1}_{3,4} \; \left(- \lambda + \nu_{y_{3,3}} \; \mu \right) + C^{k+1}_{3,3} \; \left(1 + 4 \; \lambda + \nu^{(1,0)}_{{x}_{3,3}} + \nu^{(0,1)}_{{y}_{3,3}}\right) + C^{k+1}_{2,3} \; \left(- \lambda - \nu_{x_{3,3}} \; \mu \right) + C^{k+1}_{3,2} \; \left(- \lambda - \nu_{y_{3,3}} \; \mu \right) \; = \\ \; C^{k}_{4,3} \; \left(\lambda - \nu_{x_{3,3}} \; \mu \right) + C^{k}_{3,4} \; \left(\lambda - \nu_{y_{3,3}} \; \mu \right) + C^{k}_{3,3} \; \left(1 - 4 \; \lambda - \nu^{(1,0)}_{{x}_{3,3}} - \nu^{(0,1)}_{{y}_{3,3}}\right) + C^{k}_{2,3} \; \left(\lambda + \nu_{x_{3,3}} \; \mu \right) + C^{k}_{3,2} \; \left(\lambda + \nu_{y_{3,3}} \; \mu \right) $$
Our goal was to rewrite the ecuations system, in the form:
$$ \underline{\underline{A}} \; \underline{C^{k+1}} \; = \; \underline{\underline{B}} \; \underline{C^{k}} $$
So, in that moment, we could finally put all that in a matrix form, being all $\;C^{k+1}_{i,j}\;$ the unknowns, and $\;C^{k}_{i,j}\;$ and rest of the constants, the knowns, and making zero all terms that contained i=0 or i=4 or j=0 or j=4, that are our boundary conditions. Renamed $\;\alpha \; = \; 1 + 4 \; \lambda\;$, $\;\beta \; = \; 1 - 4 \; \lambda\;$ and $\;V_{i,j}\;=\;\nu^{(1,0)}_{{x}_{i,j}} + \nu^{(0,1)}_{{y}_{i,j}}$:
$ \begin{bmatrix}
\begin{array}{ccc|ccc|ccc}
\alpha + V_{1,1} & - \lambda + \nu_{y_{1,1}} \; \mu & 0 & - \lambda + \nu_{x_{1,1}} \; \mu & 0 & 0 & 0 & 0 & 0\\
- \lambda - \nu_{y_{1,2}} \; \mu & \alpha + V_{1,2} & - \lambda + \nu_{y_{1,2}} \; \mu & 0 & - \lambda + \nu_{x_{1,2}} \; \mu & 0 & 0 & 0 & 0\\
0 & - \lambda - \nu_{y_{1,3}} \; \mu & \alpha + V_{1,3} & 0 & 0 & - \lambda + \nu_{x_{1,3}} \; \mu & 0 & 0 & 0\\
\hline\\
- \lambda - \nu_{x_{2,1}} \; \mu & 0 & 0 & \alpha + V_{2,1} & - \lambda + \nu_{y_{2,1}} \; \mu & 0 & - \lambda + \nu_{x_{2,1}} \; \mu & 0 & 0\\
0 & - \lambda - \nu_{x_{2,2}} \; \mu & 0 & - \lambda - \nu_{y_{2,2}} \; \mu & \alpha + V_{2,1} & - \lambda + \nu_{y_{2,2}} \; \mu & 0 & - \lambda + \nu_{x_{2,2}} \; \mu & 0\\
0 & 0 & - \lambda - \nu_{x_{2,3}} \; \mu & 0 & - \lambda - \nu_{y_{2,3}} \; \mu & \alpha + V_{2,3} & 0 & 0 & - \lambda + \nu_{x_{2,3}} \; \mu\\
\hline\\
0 & 0 & 0 & - \lambda - \nu_{x_{3,1}} \; \mu & 0 & 0 & \alpha + V_{3,1} & - \lambda + \nu_{y_{3,1}} \; \mu & 0\\
0 & 0 & 0 & 0 & - \lambda - \nu_{x_{3,2}} \; \mu & 0 & - \lambda - \nu_{y_{3,2}} \; \mu & \alpha + V_{3,2} & - \lambda + \nu_{y_{3,2}} \; \mu\\
0 & 0 & 0 & 0 & 0 & - \lambda - \nu_{x_{3,3}} \; \mu & 0 & - \lambda - \nu_{y_{3,3}} \; \mu & \alpha + V_{3,3}
\end{array}
\end{bmatrix}\begin{bmatrix}
C^{k+1}_{1,1}\\C^{k+1}_{1,2}\\C^{k+1}_{1,3}\\\hline\\C^{k+1}_{2,1}\\C^{k+1}_{2,2}\\C^{k+1}_{2,3}\\\hline\\C^{k+1}_{3,1}\\C^{k+1}_{3,2}\\C^{k+1}_{3,3}\end{bmatrix}$$=$
$\begin{bmatrix}
\begin{array}{ccc|ccc|ccc}
\beta - V_{1,1} & \lambda - \nu_{y_{1,1}} \; \mu & 0 & \lambda - \nu_{x_{1,1}} \; \mu & 0 & 0 & 0 & 0 & 0\\
\lambda + \nu_{y_{1,2}} \; \mu & \beta - V_{1,2} & \lambda - \nu_{y_{1,2}} \; \mu & 0 & \lambda - \nu_{x_{1,2}} \; \mu & 0 & 0 & 0 & 0\\
0 & \lambda + \nu_{y_{1,3}} \; \mu & \beta - V_{1,3} & 0 & 0 & \lambda - \nu_{x_{1,3}} \; \mu & 0 & 0 & 0\\
\hline\\
\lambda + \nu_{x_{2,1}} \; \mu & 0 & 0 & \beta - V_{2,1} & \lambda - \nu_{y_{2,1}} \; \mu & 0 & \lambda - \nu_{x_{2,1}} \; \mu & 0 & 0\\
0 & \lambda + \nu_{x_{2,2}} \; \mu & 0 & \lambda + \nu_{y_{2,2}} \; \mu & \beta - V_{2,2} & \lambda - \nu_{y_{2,2}} \; \mu & 0 & \lambda - \nu_{x_{2,2}} \; \mu & 0\\
0 & 0 & \lambda + \nu_{x_{2,3}} \; \mu & 0 & \lambda + \nu_{y_{2,3}} \; \mu & \beta - V_{2,3} & 0 & 0 & \lambda - \nu_{x_{2,3}} \; \mu\\
\hline\\
0 & 0 & 0 & \lambda + \nu_{x_{3,1}} \; \mu & 0 & 0 & \beta - V_{3,1} & \lambda - \nu_{y_{3,1}} \; \mu & 0\\
0 & 0 & 0 & 0 & \lambda + \nu_{x_{3,2}} \; \mu & 0 & \lambda + \nu_{y_{3,2}} \; \mu & \beta - V_{3,2} & \lambda - \nu_{y_{3,2}} \; \mu\\
0 & 0 & 0 & 0 & 0 & \lambda + \nu_{x_{3,3}} \; \mu & 0 & \lambda + \nu_{y_{3,3}} \; \mu & \beta -V_{3,3}
\end{array}
\end{bmatrix}\begin{bmatrix}
C^{k}_{1,1}\\C^{k}_{1,2}\\C^{k}_{1,3}\\\hline\\C^{k}_{2,1}\\C^{k}_{2,2}\\C^{k}_{2,3}\\\hline\\C^{k}_{3,1}\\C^{k}_{3,2}\\C^{k}_{3,3}\end{bmatrix}$
3) Now, we wrote the matrixes in $\;n^2\;$ matrix blocks of $\;n\;$ x $\;n\;$ dimension, and the vectors in $\;n\;$ vector blocks of $\;n\;$ dimension:
$$ \underline{\underline{A}} \; = \begin{bmatrix}
\underline{\underline{A_{1,1}}} & \underline{\underline{A_{1,2}}} & \cdots & \underline{\underline{A_{1,n}}}\\
\underline{\underline{A_{2,1}}} & \underline{\underline{A_{2,2}}} & \cdots & \underline{\underline{A_{2,n}}}\\
\vdots & \vdots & \ddots & \vdots \\
\underline{\underline{A_{n,1}}} & \underline{\underline{A_{n,2}}} & \cdots & \underline{\underline{A_{n,n}}}\\
\end{bmatrix} \;\;\;\;
\underline{C^{k+1}} \; = \begin{bmatrix}
\underline{C^{k+1}_{1}} \\
\underline{C^{k+1}_{2}} \\
\vdots \\
\underline{C^{k+1}_{n}} \\
\end{bmatrix} \;\;\;\;
\underline{\underline{B}} \; = \begin{bmatrix}
\underline{\underline{B_{1,1}}} & \underline{\underline{B_{1,2}}} & \cdots & \underline{\underline{B_{1,n}}}\\
\underline{\underline{B_{2,1}}} & \underline{\underline{B_{2,2}}} & \cdots & \underline{\underline{B_{2,n}}}\\
\vdots & \vdots & \ddots & \vdots \\
\underline{\underline{B_{n,1}}} & \underline{\underline{B_{n,2}}} & \cdots & \underline{\underline{B_{n,n}}}\\
\end{bmatrix} \;\;\;\;
\underline{C^{k}} \; = \begin{bmatrix}
\underline{C^{k}_{1}} \\
\underline{C^{k}_{2}} \\
\vdots \\
\underline{C^{k}_{n}} \\
\end{bmatrix} $$
And finally, we realised the pattern it followed, and could generalized for a $\;n\;$ x $\;n\;$ space steps:
$$ \underline{\underline{A_{w,w}}} \; = \begin{bmatrix}
1 + 4 \; \lambda + \nu^{(1,0)}_{{x}_{w,1}} + \nu^{(0,1)}_{{y}_{w,1}} & - \lambda + \nu_{y_{w,1}} \; \mu & 0 & \cdots & 0\\
- \lambda - \nu_{y_{w,2}} \; \mu & 1 + 4 \; \lambda + \nu^{(1,0)}_{{x}_{w,2}} + \nu^{(0,1)}_{{y}_{w,2}} & - \lambda + \nu_{y_{w,2}} \; \mu & \cdots & 0\\
0 & - \lambda - \nu_{y_{w,3}} \; \mu & 1 + 4 \; \lambda + \nu^{(1,0)}_{{x}_{w,3}} + \nu^{(0,1)}_{{y}_{w,3}} & \cdots & 0\\
\vdots & \vdots & \vdots & \ddots & \vdots\\
0 & 0 & 0 & \cdots & 1 + 4 \; \lambda + \nu^{(1,0)}_{{x}_{w,n}} + \nu^{(0,1)}_{{y}_{w,n}}\\
\end{bmatrix} $$
$$ \underline{\underline{A_{w,w-1}}} \; = \begin{bmatrix}
- \lambda - \nu_{x_{w,1}} \; \mu & 0 & 0 & \cdots & 0\\
0 & - \lambda - \nu_{x_{w,2}} \; \mu & 0 & \cdots & 0\\
0 & 0 & - \lambda - \nu_{x_{w,3}} \; \mu & \cdots & 0\\
\vdots & \vdots & \vdots & \ddots & \vdots\\
0 & 0 & 0 & \cdots & - \lambda - \nu_{x_{w,n}} \; \mu \\
\end{bmatrix} $$
$$ \underline{\underline{A_{w-1,w}}} \; = \begin{bmatrix}
- \lambda + \nu_{x_{w-1,1}} \; \mu & 0 & 0 & \cdots & 0\\
0 & - \lambda + \nu_{x_{w-1,2}} \; \mu & 0 & \cdots & 0\\
0 & 0 & - \lambda + \nu_{x_{w-1,3}} \; \mu & \cdots & 0\\
\vdots & \vdots & \vdots & \ddots & \vdots\\
0 & 0 & 0 & \cdots & - \lambda + \nu_{x_{w-1,n}} \; \mu \\
\end{bmatrix} $$
4) The last step focuses on isolating the unknowns vector, by doing:
$$ \underline{C^{k+1}} \; = \; \underline{\underline{A^{-1}}} \; \underline{\underline{B}} \; \underline{C^{k}} \; = \; \underline{\underline{\mathcal{M}}} \; \underline{C^{k}}$$
Where $\;\underline{\underline{\mathcal{M}}} \; = \; \underline{\underline{A^{-1}}} \; \underline{\underline{B}}$
That is the last ecuation system must be solved each time iteration. Note that, for the first iteration $\;(k\;=\;0)$, we get $\;\underline{C^{1}}\;=\;\underline{\underline{\mathcal{M}}} \; \underline{C^{0}}\;$, and this $\;\underline{C^{0}}\;$ vector, is given with the initial conditions, in our case, a Gaussian Distribution, centered in the middle of the plate.
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.