MODEL

MATLAB

Box2D

Extended Janulevicius Model

Introduction

According to professors Paul Freemont and Richard Kitney, a truthful answer to the question “What makes synthetic biology so unique?” is: “The application of engineering principles and systematic design tools involving computational modelling, modular parts and standardised measurements”[1]. This section is dedicated to the computational modelling component of synthetic biology. Although most related models aim at the prediction and understanding of the mechanisms of genetic processes, other types of models can be of great use to synthetic biology projects.

Our project, as well as structural colour itself, is characterized by order. The bacteria must be arranged in a specific, complex structure in order to produce structural colour. Therefore, the most important variables in our project are the positions of the bacteria with respect to time. The analysis of spatial and temporal coordinates is a field where mathematics and physics thrive. Therefore, in order for one to understand the mechanisms and the parameters of these spatial structures, one has to study the mathematics behind them.

As with every biological mathematical entity, the complex spatial coordination of bacteria poses a high dimensionality problem. Thus, numerous techniques which aim to reduce the dimensions of this problem have been developed[2].

Continuous models do not take into account properties of each individual cell, but of the average quantity of them at a particular point in space. They have relatively few degrees of freedom and they tend to be computationally efficient. However, they are unable to evaluate the specific coordinates of each cell, and are therefore unusable in predicting complex structures[2].

Cellular Potts Model/Glazier-Graner-Hogeweg (GGH) Model[3], which is a stochastic on-lattice model. The cells occupy numerous lattice sites and the dimensions and flexibility of each cell can be computed using a Monte-Carlo simulation approach. The problem with this type of models is their high computational expense and mechanical inaccuracy.[2]

Node-based Monte-Carlo methods, which lack direct correlation with experimental measures.[2]

In addition to these classical approaches, numerous other specialized models have been proposed, such as the one proposed in this essay.

It has been observed that the most important factor in the development of the ordered structure that Flavobacteriia exhibit is their special manner of movement, namely, gliding motility[4]. A bacterium that exhibits similar type of movement is the proteobacterium Myxococcus xanthus, which has been extensively modelled. After examining and evaluating various gliding motility models, we decided to develop a model based on the one proposed by Janulevicius et al.[5] (hence in this study referred to as “Janulevicius Model”). This model was chosen as the basis of this analysis due to its flexibility to changes, its ability to predict the exact positions of each bacterium, its computational efficiency, the intuitive way it is built, its suggested potential to be extended into a model for Flavobacteriia[6] and the fact that its parameters, being purely biophysical, can be experimentally measured. The Janulevicius Model is a mass-spring model, and it is specifically made to predict the movement of M. xanthus. The main difference between the movement of M. xanthus and F. johnsoniae is that the latter also exhibit a rotational movement in addition to their linear one[7]. Therefore, the following model we developed is an extended version of the Janulevicius Model, in which an additional rotational engine force has been applied

The Model

A system of a number of bacteria, M, is assumed. Each bacterium, of length L and width W, is modelled as an ordered array of N particles, which are bound together via linear and angular springs. Each particle i = 1,...,N of every bacterium j = 1, … ,N is characterized by its position \(\vec{r}_{i}\), and velocity \(\vec{u}=\dot{\vec{r}}_{i}\), and by the forces that are acted upon it,\(\vec{F}_{i}\) . The force of each particle, \(\vec{F}_{i}\), is the sum of the following components:

Linear Spring Forces

They are the forces that resist elongation or shortening of the bacterium. Every two adjacent particles of a bacterium are considered to be bound by a linear spring, which is defined by the vector:\[\vec{l}_{i j} \equiv \vec{r}_{i+1} - \vec{r}_{i} \quad \forall i=1, \ldots, N-1, (1)\]
Each linear spring of a bacterium acts to each adjoining particle the force: \[F^{\vec{l}^{l, i}}_{i j}=-k^{l} \cdot\left(l_{i}-l_{0}\right) \cdot\left(\frac{\vec{l}_{i}}{l_{i}}\right) \quad \forall i=1, \ldots, N-1 ,(2) \] \[F^{\vec{l}^{l, i}}_{(i+1) j}=+k^{l} \cdot\left(l_{i}-l_{0}\right) \cdot\left(\frac{\vec{l}_{i}}{l_{i}}\right) \quad \forall i=1, \ldots, N-1 , (2)\]
\(k^{l}\) is the stiffness of the spring.Its value increases as the resistance to elongation and shortening increases.

Angular Spring Forces

They are the forces that resist bending of the bacterium. Every three adjacent particles of a bacterium are considered to be bound by an angular spring, which is defined by the angle between the vectors \(\vec{l}_{i j} \text { and } \vec{l}_{(i+1) j}, \alpha_{i j}, \forall i=1, \ldots, N-2\). Each angular spring produces two torques: \[\vec{\tau}_{i j}^{a, i}=k^{a} \cdot \alpha_{i j} \cdot\left(\frac{\vec{m}_{i j}}{m_{i j}}\right) \quad \forall i=1, \ldots, N-2 ,(3)\] \[ \vec{\tau}_{(i+2) j}^{a, i}= -k^{a} \cdot \alpha_{i j} \cdot\left(\frac{\vec{m}_{i j}}{m_{i j}}\right) \quad \forall i=1, \ldots, N-2,(3) \] where \(\vec{m}_{i j} \equiv \vec{l}_{i j} \times \vec{l}_{(i+1) j} \nonumber\). These torques produce the forces: \[\overrightarrow{F^{a}}_{i j}^{\alpha, i}=\frac{1}{l_{i j}^{2}} \cdot\left(\vec{l}_{i j} \times \vec{t}_{i j}^{a, i}\right),(4)\] \[\overrightarrow{F a}_{(i+2) j}^{\alpha, i}=-\frac{1}{l_{(i+1) j}^{2}} \cdot\left(\vec{l}_{(i+1) j} \times \vec{\tau}_{(i+2) j}^{a, i}\right) \quad \forall i=1, \ldots, N-2,(4) \] \[ \overrightarrow{F^{a}}_{(i+1) j}^{\alpha, i}=-\left(\overline{F^{a}}_{i j}^{\alpha, i}+\overrightarrow{F^{a}}_{(i+2) j}^{\alpha, i}\right),(4)\]
\(k^{a}\) is the stiffness of the spring.Its value increases as the resistance to bending increases and the bacterium becomes more rigid.

Linear Engine Forces

Gliding motility [also known as adventurous (A) motility], is defined as the active surface translocation along the long cell axis, without the aid of flagella or pili, and it is responsible for individual cell movement[6]. The bacterium glides at constant velocity, \(u^{b}\). The mechanisms of gliding are not well understood. Gliding in different bacteria involves a broad spectrum of underlying mechanisms[8]. In F. johnsoniae, the function of the gliding motor is monitored through the motion of SprB and RemA[8]. It should be noted that although the biological mechanisms of the gliding motor differ between the two species, their movement is comparable. Therefore, the linear engine of M. xanthus may be extrapolated to F. johnsoniae. In this model, gliding motility is stimulated by a linear engine. The force of the engine is considered to be evenly distributed among the particles of each bacterium. The magnitude of the force,\(F^{e}\) , can be measured experimentally. Thus, the force acted on each particle of each bacterium is equal to: \[\vec{F}_{ij}^e=\left(\frac{F^e}{N}\right)\cdot{\hat{t}}_{ij} \quad \forall i=1, \ldots, N ,(5)\] \[\text{where } {\hat{t}}_{ij}\equiv-k^e\cdot\frac{\vec{l}_{(i-1) j}+\vec{l}_{ij}}{\lVert{\vec{l}_{(i-1) j} + \vec{l}_ij}\rVert} \quad \forall i=2, \ldots, N - 1 \] \[{\hat{t}}_{1j}\equiv -k^e\cdot \frac{(\vec{l}_{1j})}{\lVert\vec{l}_{1j}\rVert} \quad \forall i=2, \ldots, N - 1 \] \[{\hat{t}}_{Nj}\equiv -k^e \cdot \frac{(\vec{l}_{(N-1) j})}{\lVert\vec{l}_{(N-1) j}\rVert} \quad \forall i=2, \ldots, N - 1\]
\(k^{e}\) is a constant that can take the values -1 and +1 and it gives the direction of the engine. The direction of a bacterium is reversed if \(t-t_{LR} > T_R\), where t is the current time of the simulation, \(t_{LR}\) is the time of the last reversal of the bacterium, and \(T_R\) is time interval until the next reversal.

Rotational Engine Forces

As stated before, one of the main differences between the movement of M. xanthus and F. johnsoniae is that the latter also rotates with one of its ends fixed. Therefore, a rotational force has to be applied to the model. Solving Newton’s equations with the constraint that every particle of the same bacterium has the same angular velocity, yields that for every bacterium j, the following variable is constant: \[c_{j} \equiv \frac{\left\lVert\vec{F}_{i j}^{-r}\right\rVert \cdot \sin \left(\theta_{i j}\right)}{\left\lVert\vec{r}_{i j}\right\rVert},(6)\]
where \(θ_{ij}\) is the angle between the vector of the stationary particle of the bacterium, and the particle i. Each value \(c_j\) can be calculated by the formula: \[c_j=\frac{\tau}{\sum_{(i\neq stationary\ end)}r_{ij}^2},(7)\] where 𝜏 is the torque of the rotational engine.Therefore, \(\lVert\vec{F}_{ij}^{a}\rVert\) can be easily calculated from the equation (6). The direction of the force is given as: \[{\vec{F}}_{ij}^r=k^r\cdot\lVert\vec{F}_{ij}^{a}\rVert\cdot\hat{n}_{ij} \quad \forall i\neq stationary\ end\,(8)\]
where \(\hat{n}_{i j} \perp \hat{t}_{i j}\).\(k^{r}\) is a constant that can take the values -1 and +1 and it gives the direction of the engine.

Collision Forces

Collision Detection

Different bacteria (or different parts of the same bacterium) can overlap, resulting in collision. For the detection of the collision, each bacterium is considered to be an ordered array of line segments, defined parametrically as \(\vec{Q}_{i j}(P) \equiv \vec{r}_{i j}+P \cdot\left(\vec{r}_{(i+1) j}-\vec{r}_{i j}\right)\), where \(P \in[0,1]\). A collision is considered to be taking place if the magnitude of the distance between two non-adjacent line segments and becomes smaller than the critical value W. The algorithm used to find the variables \(\vec{d}, P_{1} \text { and } P_{2}\) and can be found in the book “Real-Time Collision Detection” by C. Ericson[9].

Collision Resolution

When a collision takes place between the line segment \(Q_{ij}\) and \(Q_{kl}\), the following forces are applied to the particles: \[{\vec{F}^c}_{ij}=-\left(1+P_1\right)\cdot\left[k^c\cdot\left(d-W\right)\cdot\left(\frac{\vec{d}}{d}\right)\right] \quad \forall i=1, \ldots, N - 2,(9) \]
\[{\vec{F}^c}_{(i+1) j}=-P_1 \cdot\left[k^c\cdot\left(d-W\right)\cdot\left(\frac{\vec{d}}{d}\right)\right] \quad \forall i=1, \ldots, N - 2,(9) \]
\[{\vec{F}^c}_{kl}=-\left(1+P_2\right)\cdot\left[k^c\cdot\left(d-W\right)\cdot\left(\frac{\vec{d}}{d}\right)\right] \quad \forall i=1, \ldots, N - 2 ,(9)\]
\[{\vec{F}^c}_{(k+1) l}=P_2 \cdot\left[k^c\cdot\left(d-W\right)\cdot\left(\frac{\vec{d}}{d}\right)\right] \quad \forall i=1, \ldots, N - 2,(9)\]
\(k_c \) is the stiffness of the collisions. Its value is chosen so as to prevent overlapping of bacteria. In addition, in order to prevent excessive bending of angular springs, the following forces are introduced: \[ {\vec{F}^c}_{ij}=-k^c\cdot\left(d-W\right)\cdot\left(\frac{\vec{d}}{d}\right) \quad \text { If } d \equiv \vec{r}_{i j}-\vec{r}_{(i+1) j} \leq W,(10) \] \[ {\vec{F}^c}_{(i+1) j}=-k^c\cdot\left(d-W\right)\cdot\left(\frac{\vec{d}}{d}\right) \quad \text { If } d \equiv \vec{r}_{i j}-\vec{r}_{(i+1) j} \leq W,(10)\]

Equations of Motion

It is assumed that inertia effects are negligible and that the velocity of each particle is equal to its terminal velocity at every time[5]. This yields the following system of ordinary differential equations (ODEs) that is applied for every particle of every bacterium:
\[ \frac{d \vec{r}_{i j}(t)}{d t} \equiv \vec{u}_{i j}=\left(\frac{1}{\zeta^{t}}\right) \cdot\left(\hat{t}_{i j} \cdot \vec{F}_{i j}\right) \vec{F}_{i j} +\left(\frac{1}{\zeta^{n}}\right) \cdot\left(\hat{n}_{i j} \cdot \vec{F}_{i j}\right) \vec{F}_{i j}, (11) \] \[\text{where } \zeta^t = \frac{(\frac{F^e}{N})}{u_b} \text{, } \zeta^n=2\zeta^t \nonumber \]
In this case, the bacteria move in a two-dimensional space. Therefore, the system (11) consists of \( 2 \times N \times M \) ODEs.

The system (11) is solved numerically using the variable order method. As stated in the Contribution section, a custom MATLAB function has been developed that computes the values of the system (11) for each time step. You can find that file here.

The parameters, along with their sources, are listed in Table 1.

Parameters	Source
\( L = 3.98 \mu m \)	4
\( W = 298 nm \)	4
N = 8
\( k^l = 10^{-2} \)	Assumed
\( k^a = 10^{-15} Nm \)	Assumed
\( k^c = 5.5 \times 10^{-4} NM \)	Assumed
\( F^e = 100 pN \)	Assumed
\( \tau = 1000 pN * mm \)	7
\( k^{r}=\left\{\begin{array}{l} +1,92 \% \text { probability} \\ -1,8 \% \text { probability } \end{array}\right. \)	7
\( u_b = 3 \mu m / s \)	7

Analysis

From visual observation of the simulation plots, it is clear that the model is successful in qualitatively predicting the movement of the bacteria. Moreover, the simulator that has been developed is computationally stable in most cases of initial values.

It is noted that due to the lack of experiments, it was impossible to get accurate values for the parameters. Therefore, some of them were found from the literature, while others were assumed (Table 1) based on optical observation of movies found in scientific papers. However, during the literature search, it was found that most of these parameters change significantly for each different strain of F. johnsoniae. As a result, it is not possible to get the exact prediction of the movement of the strains that are examined in this project. That being said, the team is going to extend the project for the iGEM 2021 competition and then, hopefully under the current COVID-19-limiting circumstances, the proper values will be known experimentally from the wet lab.

Nonetheless, some very important qualitative results can be obtained. Two cases are simulated. In the first case (Movie 1), the torque of the bacteria is zero. In this case, as expected, the result is similar to the one obtained by the original Janulevicius model (but with different values of the parameters). On the other case (Movie 2), the torque value presented in Table 1 is applied. As it can been seen from the Movie 2, the torque changes significantly the direction of the bacteria after their collision. Therefore, the implementation of the extra rotational force in the Janulevicius model alters significantly the arrangement of the bacteria. This alteration can be an important factor in the construction of the complex arrangement of the Flavobacteriia, that ultimately produces structural colour. Similar differences can be observed for different initial conditions and number of bacteria in the system.

Thoughts about the Future

The Janulevicius model has been extended by its developer by making more assumptions for the biophysics of the M. xanthus bacteria [2]. It is important to examine whether these [7] assumptions are true for F. johnsoniae before applying them to this model. Furthermore, additional kinds of movement have been reported for F. johnsoniae [10].It is apparent that the model can become more accurate via the implementation of more relevant forces. Nevertheless, the selection and development of this model and its simulator had been a time-consuming process. Moreover, in order to properly introduce these forces into the model, detailed experimental data are needed. Consequently, one of the main goals of the iGEM Athens team is to expand and further improve the model and the simulator in the next phase of the iGEM competition.

Gliding Motility Simulation

In order to simulate the gliding motility mechanism that Flavobacteria have, we created a simulation based on Newtonian interactions using Box2D, a 2D physics engine for games. To study the mechanical cellular interactions, we examined the case of head-to-side collision between two Flavobacteriia moving in perpendicular directions and then we scaled it up to a large population of Flavobacteria moving in different random directions. The simulation was based on [6] and [5] models for Myxoccocus xanthus and is in progress. We investigated the case in which Flavobacterium johnsoniae are represented as rigid bodies and as flexible bodies, respectively. Regarding the case of the flexible bacteria, two types of models were constructed, the Viscous Coupling Model (VCM) and the Elastic Coupling Model (ECM). To differentiate them we assumed strong attachments between the cells and the substrate in the ECM. This assumption leads to the primary cell (the one whose side is hit) maintaining its direction after the collision and the secondary cell aligning with the former. Respectively, in the VCM model, after the collision, both cells changed directions and aligned to a new common one.

Each cell is modelled as a set of connected circles (nodes) and rectangular spacers. On the one hand, supposing that the bacteria are rigid bodies, each cell is modelled as 2 circular nodes connected by a rectangular spacer. On the other hand, supposing that the bacteria are flexible bodies, but still rather stiff, each cell is modelled as 7 circles and 6 rectangles in the Box2D world, connected by rotational joints and distance joints. In both cases, we set the radii of the circles and both the width and length of the rectangles to appropriate values in an attempt to replicate the total length:width ratio of F. johnsoniae [11], since we can’t work with the actual dimensions that are in the micron scale, via the Box2D environment.

In the case of the bacteria being rigid bodies, the construction of each cell is much simpler and it only uses 2 rotational joints connecting each circle to the rectangular spacer. For the rotational joints we used Box2D’s built in Rotational Joints and set the angle boundary to [-π/8, π/8] to keep the cell from bending too much. The rotational joint in this case only acts as a means to stick all the elements together.

For the simulation of the bacteria being flexible bodies, we used one rotational joint between the circles and the rectangular spacers between them, creating 12 rotational joints.

Figure 2: Shape of modeled Flavobacteriia as flexible bodies (top) and as rigid bodies (bottom)
Rotational Joints are marked as red dots.

We also used Box2D’s built in Distance Joints and set its parameters to proper values in order for it to act like a very stiff spring and its length to 2*h, where h is the height of every rectangular spacer, to connect every circular node (i) to the one after its next one (i+2).The function of the latter is to control the bending of the cells by keeping the distance between nodes indexed i and i+2 fixed, and making sure that the cell will go back to its original shape after every collision. Additionally, we added a counteracting force to every node, in order to restore any displacement caused by collisions. This force simulates a linear spring between the node (i) and the middle point of the line segment that connects the node to its left (i-1) and the one to its right (i+1).

Figure 3: Straightening force applied to displaced nodes of flexible Flavobacteriia.

In the Elastic Coupling Model we also implemented the (strong) node-substrate interaction force that is applied to every circular node when it is displaced and restores it to its original position. Without the assumption of strong node-substrate interaction, the implemented simulation refers to the Viscous Coupling Model.

Figure 4: Node-substrate interaction force applied on displaced nodes of flexible Flavobacteriia.

In both the proposed mechanisms of gliding motility, propulsive force is generated by the cell’s distributed “engine”. We assumed that any forces on the movement axis along with the cell motility force result in a constant linear velocity equal to the actual mean velocity of F. johnsoniae (4μm/s [11]). In Box2D, assuming flexible cells, we apply constant velocity vectors (initial condition) to all nodes except the first and the last ones thus simulating the above. Assuming rigid cells, we apply force to the middle of the rectangular spacer.

Figure 5: Velocity initial condition applied to flexible (top) and rigid (bottom) cells.

As for the collisions, Box2D handles them by default, not letting bodies overlap with each other and applying the appropriate collision resolving forces. Box2D also solves the equations of motion for all bodies and outputs the position and velocity at each time step.

Scaling up from 2 cells to 500 cells we captured these videos for both the flexible and the rigid cases, for the Viscous Coupling Model. All the parameters were chosen both form literature and via trial and error process to match our estimations for the time being.

Figure 6: Top Left: Simulation of 2 flexible cells using the Elastic Coupling Model (top left)
Top Right: Simulation of 2 flexible cells using the Viscous Coupling Model (top right)
Bottom Left: Sped up Simulation of 500 flexible cells using the Viscous Coupling Model (bottom left)
Bottom Right: Sped up Simulation of 500 rigid cells using the Viscous Coupling Model (bottom right).

The improvement of the described simulation is still in progress. It will be fine-tuned to match F. johnsoniae, that we plan on using in our project, once we acquire enough experimental data. The parameters will be fixed to better match reality and more elements are going to be added (like the spontaneous turning of the cells [7]). Also, one of the two cell models (flexible/rigid) will finally be eliminated based on our lab data and the simulation will be further optimised.

Team:Athens/Model