Team:Costa Rica/Model

Mathematical modelling is a tool used to predict the behavior of a problem without needing endless amounts of experimentation. In this section, we present the equations required to accurately describe our problem.

Index

Introduction
Predator-Prey Model
Conclusions
Code
References

Introduction

Exponential & Logistic Growth

It is known that living organisms, and specially microorganisms tend to reproduce ex- ponentially, under the assumption of unlimited resources and space, and in the abscense of predators or catastrophic events. In the real world, food is inevitably going to run out as the number of individuals grows, this causes a slowdown in the exponential growth. This is the main and most basic principle that models the growth and behavior of a single species population, under the realistic assumption of limited resources over time.

Mathematically, exponential growth can be expressed in the following way: we have a number of individuals $x$ which is going to change over time $t$. In this case, it is known that a higher number of individuals means a higher rate of population growth. The growth rate of our population is going to be written as $x'$ (the derivative of $x$ over time). It is now possible to model this type of growth via a simple equation:

(1.1)

\begin{equation}\label{eq:1} x' = Kx \end{equation}

Where $K$ is some proportionality constant. This says that the population's growth rate is directly proportional to the number of living members. Equation (1.1) is a differential equation, since it relates a certain quantity and rate of change over time (its derivative). The solution of this equation is going to predict the exact number of individuals as a function of time:

(1.2)

\begin{equation} x(t) = Ce^{Kt} \end{equation}

Where $C$ is a constant that arises in the process of solving the equation. On specific circumstances, like when the starting population number is given, $C$ can be easily calculated.

Figure 1 shows the behavior of a theoretical population modeled by equation (1.1)

Figure 1. Exponential Model.
Source: MATLAB generated.

In this theoretical population, an initial number of only 5 individuals grows to a number of 20 million in 15 time units. If more time units were to be simulated, the population number will tend to infinity very quickly. It is reasonable to conclude that no real world population is going to develop in this fashion. Whether it is by natural causes, like lack of resources or space (or competition with other organisms for these), or by external causes, like random catastophic events, the number of individuals $x$ cannot grow uncontrollably. Therefore, it is necessary to improve our modeling tools.

The main reason the exponential model is presented is that, although it is not realistic by any means, it is the first step in formulating better and more accurate mathematical models

We are now going to take into account the simple fact that an environment cannot have infinite resources, since by the previous discussion, a population simply cannot grow towards infinity. The logistic model is now introduced, which considers the situation in where there cannot be too many individuals coexisting on a limited amount of resources or space.

(1.3)

\begin{equation} x'(t) = x(1-x) \end{equation}

In its simplest form, $x$ represents once again the number of individuals of our population and, as before, the more individuals, the higher the growth rate ($x'$). Nevertheless, the term $(1-x)$ works as a limiting factor, it inhibits the growth rate by a certain amount, which is propotional to the number of individuals. An intuitive way of understanding this situation is the following: as $x$ gets bigger, it is going to approach a ''maximum capacity" which is in some way, a $100 \%$ of sustainable population with our given resources. This way, the term $(1-x)$ means how close is $x$ to not being able to sustain itself with the resources that the environment can provide. When $x$ is too close to $1$ (which in this case is called the carrying capacity), the $(1-x)$ term is going to be very close to $0$, effectively making $x'$ very small. This is translated as: ''Too many individuals will slow down growth''. If equation ($1.3$) is solved, the function obtained will be:

$$x(t) = \frac{e^t}{e^t + C}$$

Where $C$ can again be calculated if the initial conditions of the problem are known.

Figure 2. Logistic Model.
Source: MATLAB generated.

When $x(t)$ is graphed, a classic sigmoid (''S'' shaped) curve is obtained. These curves are usually encountered when there is a factor that inhibits growth when $x$ is high (not only encountered in biology). For this population, the starting population represented $5 \% $ of the carrying capacity (it is necessary to work in percentages for this model). It can be seen that, innitialy, since there are too few individuals, growth is slow. As $x$ increases, at about $1/3$ of the total capacity, there is a demographic boom, since there are now enough individuals and resources to grow succesfully. As the population number reaches $100 \%$ of the system capacity, it can be seen in the graph that the slope of the function (or the growth rate) slows drastically until it reaches some kind of equilibrium. A population that has reached a sustainability limit will not grow past it (according to this model at least).

For the purposes of our study, a more elaborated logistic model is presented, where the carrying capacity is no longer $1$ but some number $K$. This way, the term $x/K$ still means what fraction of $K$ does $x$ represent. Moreover, the reproduction rate $r$ is now introduced, which controls more closely the exact proportion between $x'$ and $x$. Both $K$ and $r$ should be calculated empirically (or taken from literature), since they depend completely on the modelled species, and its environment. Our base equation is now written as

(1.4)

\begin{equation} x' = rx\left(1-\frac{x}{K} \right) \end{equation}

Differential Predator-Prey Model

Let $x(t)$ and $y(t)$ be the concentrations of C. difficile and L. casei, respectively, as functions of time. Our main goal is to eradicate $x(t)$ or at least to lower it to a safe number (in order to eliminate or avoid infection). Without specifying initial conditions, we propose the following model.

(2.1)

\begin{equation} \begin{cases} x'(t) = r_1 x \left(1- \cfrac{x}{K}\right) - \eta \cfrac{L(x)xy^2}{x^2+y^2+A} \\ y'(t) = r_2y\left(1-\cfrac{y}{K}\right) \end{cases} \end{equation}

Where

$$L(x) = \frac{Px^{\alpha}}{1+Px^{\alpha}}$$

And the constants are:

•$r_1,r_2$: Growth rates of C. difficile and L. casei, respectively. These are taken as $r_1=0.6416$ h$^{-1}$ (Thanissery et. al 2017), and $r_2=0.65$ h$^{-1}$ (Avinash, 2018), these values are standard for our lab procedures, although they can be calibrated experimentally, their order of magnitude wouldn't change significantly.

•$K$: Carrying capacity. It is fixed as a high value ($K=200$ g/L), since we want to model the L. casei versus C. difficile interaction in a rich environment, such as an intestine. We are assuming that the amount of competition for resources is not significant.

•$A$ is a mass constant, that corresponds to other mass present in the system other than biomass. For our model, we took a standard $A=20$ g/L. representing the initial concentration of glucose in our experiments.

•$P$ and $\alpha$ are Hill's equation parameters (Goutelle, 2018), and are generally taken as $P>>0$ and $\alpha =2$. In our model, $P=100$ for stability purposes.

•$\eta$ is a dimensionless, efficiency constant. It will be varied in addition to initial values to produce many possible scenarios.

We can observe, in our model (equation $(2.1)$), that the main term in each equation is the modified logistic growth term. These model each species' growth independently, as if they were not interacting at all. The last term in the fist equation of $(2.1)$, is the interaction term. This term describes how L. casei attacks C. difficile. As we know, this interaction is rather indirect, because we need to quantify the concentration of Lysine, which is produced by L. casei, in response to C. difficile's quorum sensing. In collaboration with TAS Taipei team, and according to Goutelle (2018), Hill's equation is a good tool to model Lysine concentration as a function of $x$. We multiply $L(x)$ by $xy^2$ to make it sensitive to both C. difficile and L. casei concentration, and finally divide by $(x^2 + y^2 +A)$ to normalize.

In simple words, our interaction function behaves like this: If there is too little $x$, no quorum sensing takes place, and no interaction should occur. If there is a lot of $x$, and little $y$, quorum sensing would take place, but not enough Lysine would be produced to kill C. difficile, so its concentration would not go down. Finally, if there is enough of each, both quorum sensing and Lysine production take place, effectively (depending on $\eta$) killing C. difficile population.

Remarks: Some models include a linear death rate term in each equation, but the team decided not to include it, since in literature, most of these models simulate large amounts of time, in which several life cycles of each organism occur. This does not apply to our problem, and we are more interested in modelling the effectiveness of our cure. Moreover, L. casei's life cycle is of no interest, since if it were very short, patients can always take more of the bioengineered probiotic.

Results and Scenarios

It is very likely that our ODE system cannot be solved explicitly ( or ''by hand") because of its non-linear nature. However, with the aid of software, it is possible to approximate its solution by numerical methods, and assuming stability, to obtain an accurate approximation of our model.

We now present possible scenarios, varying initial conditions and the efficiency parameter $\eta$. We emphasize that every scenario presented is qualitative, that is, what should be observed is the general behavior and shape of the curves generated, not their exact values. In order to achieve precise, quantitative results, large amounts of experimentation and calibration are needed. In every scenario, a time of $50$ h is simulated.

Nine scenarios are presented, each with three values of efficiency ($\eta$), and three initial conditions: $x\ll y$ (high initial L. Casei), $x\gg y$ (high initial C. difficile) $x=y$ (similar initial conditions).

Low Efficiency Simulations.

Setting $\eta=0.2$ ($20 \%$ efficiency), we present the possible scenarios.

Figure 3. Low efficiency, equal starting conditions ($1:1$ ratio).
Source: MATLAB generated.

Assuming lysine is not very effective, it is observed in figure 3 that, when the initial concentrations are equal, both populations will simply grow towards equilibrium, which is the maximum sustainable population $K$ allows. It should be noted that, even at low efficiency values, the interaction term is effective in inhibiting the equilibrium value for $x(t)$. If this term is omitted, both curves would be exactly the same (as in figure $2$). It is also expected that $y(t)$ will always tend to an equilibrium value, since no interaction is modeled for $y(t)$. Once again this is of little interest since if $y$ decays, the patient can take more medication.

Figure 4. Low efficiency, high initial C. difficile concentration ($20:1$ ratio).
Source: MATLAB generated.

Figure 4 describes what would happen if too much C. difficile is present at the start of the simulation. Once again $y(t)$ balances out. On the other hand, $x(t)$ has a small growing period, before being inhibited and achieving a lower than $K$ value of equilibrium.

Figure 5. Low efficiency, high initial L. casei concentration ($1:20$ ratio).
Source: MATLAB generated.

The last simulation for low efficiency lysine (fig. 5) presents a slow, inhibited growing C. difficile, that ultimately finds its equilibrium. In all three cases, the lysine only seems to inhibit the growth of $x$, and fails to avoid or eliminate the infection.

Moderate Efficiency Simulations.

Setting $\eta=0.5$ ($50 \%$ efficiency), we present the possible scenarios.

Figure 6. Moderate efficiency, equal starting conditions ($1:1$ ratio).
Source: MATLAB generated.

With equal starting concentrations, figure 6 presents a very restricted curve for $x(t)$. The equilibrium value achieved is much lower than the environment's carrying capacity, and it is reasonable to assume that any possible infection was avoided. It should be noted that in this scenario, C. difficile appears to grow in the first hours, but after a short time the lysine takes effect, succesfully limiting its growth.

Figure 7. Moderate efficiency, high initial C. difficile concentration ($20:1$ ratio).
Source: MATLAB generated.

As expected in typical logistic growth, when the population number is near half of the carrying capacity, it should grow about exponentially, which is exactly what occurs in figure 7. However, as the small concentration of L. casei also blooms, they start producing enough Lysine that quickly kills C. difficile individuals, effectively diminishing the value of $x(t)$ until it reaches a safe, infection-free equilibrium.

Figure 8. Moderate efficiency, high initial L. casei concentration ($1:20$ ratio).
Source: MATLAB generated.

Figure 8 simulates the last moderate efficiency scenario, in which there is very little C. difficile to start with. It is observed that $x(t)$ takes a very long time (almost the entire simulation) to reach its equilibrium value, which is once again, very low compared to the total capacity of the environment. The growth of the infection is greatly limited by the action of the bioengineered L. casei.

High Efficiency Simulations.

Finally, setting $\eta=0.9$ ($90 \%$ efficiency), we present the possible scenarios at once.

Figure 9. High efficiency, equal starting conditions ($1:1$ ratio).
Source: MATLAB generated.

Figure 10. High efficiency, high initial C. difficile concentration ($20:1$ ratio).
Source: MATLAB generated.

Figure 11. High efficiency, high initial L. casei concentration ($1:20$ ratio).
Source: MATLAB generated.

All three scenarios are analogous to their moderate efficiency counterparts. If the starting conditions are similar (figure 9), $x(t)$ will not be even able to grow at all, and it decreases rapidly until a very low equilibrium point (almost zero). On the other hand, when the initial value of $x(t)$ is high (figure 10), there is again a short-lived bloom of C. difficile, followed by a blook of L. casei and a sharp decrease in C. difficile population, and stabilizing again at almost zero concentration. Finally, when $x(t)$ has a low starting value compared to $y(t)$ (figure 11), its growth is severly restricted, to the point of being almost constant (next to 0), this being the most extreme scenario. It is important to note that, in figure 11, due to the extreme nature of the initial conditions and efficiency value, the last simulation does not capture correctly the quorum sensing, lysine production phenomena, as it shouldn't occur in the first place. Nevertheless, we consider this case as too extreme, and either way, no infection is present in this scenario.

Conclusions

In addition of the scenarios presented, many other possibilities were tested in order to study the system's stability. No extreme solutions were found. More specifically, no infinite, negative or wildly oscilating solutions were obtained (starting with realistic initial values), this contributes to the model credibility. As for the scenarios that were presented, at low efficiency, the lysine was not effective enough as to eradicate C. difficile individuals, although in some way, it inhibited their growth. However, in medium and high efficiency simulations, the infection was eliminated (or avoided) in both cases, the only difference being the equilibrium value for $x(t)$, which as expected, turned out to be much lower in the high efficiency scenarios.

Code

All simulations were run with the aid of MATLAB software. The following code was used for testing and calibration.

                                    

                                           %%Diff-easy. iGEM Costa Rica 2019%%

%Numerical approximation of the differential model proposed to predict
%the C. difficile vs L. casei problem
%x: C. dificile population
%y: L. casei population


%%Theorical constants%%
r1=0.6416;  %L. casei growth rate at 20g/L glucose concentration
r2=0.65;  %C. difficile growth rate (630 strain)
K=200; %Carrying capacity. Assumed to be high

eta=0.9; %Lysine efficiency

%%EMPIRICAL CONSTANTS %% 
A=20; %Mass constant (glucose level)
alpha=2; %Hill exponent 
P=100; %Hill proportionality constant

%Remark: If P>1, its value does not affect the solution considerably
%but if it is close to 0, the system may become unstable.




%%INITIAL CONDITIONS%%
t=[0,50];  %50 hour simulation time
S0=[5,100];  %S0 is the initial condition vector,
%[initial C. difficile, initial L. casei]


%The system is initialized
f = @(t,s)[r1*s(1)*(1-s(1)/K) - (eta* s(1)*(P*s(1)^alpha)/(1+P*s(1)^alpha)*s(2)^2)/(s(1)^2 + s(2)^2+ A)
    r2*s(2) * (1-s(2)/(K))];
   
[t,s] = ode45(f,t,S0);  %Numerical solution approximated using ode45 built-in command

%Graphing options
figure 
plot(t,s(:,1),'r','Linewidth',1.5)

hold on
plot(t,s(:,2),'g','Linewidth',2)

grid on
ylim([-10,250]);
xlabel('$t$','Interpreter','Latex')
ylabel('$(x,y)$','Interpreter', 'Latex')
dim = [.75 .2 .1 .2];
title('\textit{L. casei} vs. \textit{C. difficile}','Interpreter','Latex')
str = {'$\eta = 0.9$'};
legend({'$x(t)$','$y(t)$'},'Location','east','FontSize',22,'Interpreter','Latex')
annotation('textbox',dim,'String',str,'FitBoxToText','on','FontSize',20,'Interpreter','Latex');
set(gca,'FontName','Arial','FontSize',22,'FontWeight','Bold')

References

T. Avinash, Optimization of process parameters and estimation of kinetic parameters for lactic acidproduction by lactobacillus casei mtcc 1423, Biomass Conversion and Biorefinery, (2018).

S. Goutelle, The Hill equation: a review of its capabilities in pharmacological modelling, Fundamental& Clinical Pharmacology 22 (2008) 633-648, (2008).

R. Thanissery et al,Inhibition of spore germination, growth, and toxin activity of clinically relevant C. difficile strains by gut microbiota derived secondary bile acids, Anaerobe, (2017).