Discover the second part of our model, which focus on the BacProtacs part, on the Docking page!
To support laboratory tests, our team carried out a mathematical modeling study. We aimed to provide predictions on the behavior of our tool while exploring various possibilities and questions other than through manipulations - mathematical theory being far less costly!
Our model will focus on the most experimentally advanced part of our project, the CRISPR DCas9 part.
Note that the models we will present were entirely designed by our team!
SuperBugBuster aims to make bacteria lose their resistant trait to limit the spread of antibiotic resistance, the silent pandemic of our modern world. SuperBugBuster is based on two experimental approaches, one using dCas9 CRISPR technology. We use donor bacteria with a genetically modified plasmid to cause the bacteria in the environment to lose their resistance. This plasmid tool creates mutations in the gene coding for antibiotic resistance, rendering it non-functional. It is transmitted from a donor bacterium to a recipient bacterium by conjugation. For further details, go to the Description and Experiments pages.
As a reminder, the upcoming models are focused on the CRISPR dCas9 portion of our experiments. Our objective here is to comprehend how populations of sensitive and resistant bacteria will evolve in the presence of our SuperBugBuster tool. Therefore, we are turning to population dynamics models. Population dynamics mathematics are employed not only to model population growth but also to capture various interactions that can occur between them. In this case, our five populations are as follows:
\( R(t) \), population of resistant bacteria at time \(t\),
\(S(t)\), population of sensitive bacteria at time \(t\),
\(Rc(t)\), population of resistant bacteria containing the CRISPR tool plasmid at time \(t\),
\(Sc(t)\), population of sensitive bacteria that are recipients containing the CRISPR tool plasmid at time \(t\),
\(C(t)\), population of donor bacteria containing the CRISPR tool plasmid at time \(t\).
The difference between the populations \( \frac{Rc}{R} \) and \( \frac{Sc}{S} \) lies in the possession of the CRISPR tool plasmid. The population \( C \) is the donor population of the CRISPR tool plasmid.
Finally, we introduce a sixth variable: \( A(t) \), for anhydrotetracycline. This variable, our inducer, is somewhat unique, and we will revisit it in the model explanations.
Our model is initially formulated as a set of 6 differential equations. We present the final form of our model here, which is then explained further on this page.
This model serves to demonstrate the effectiveness of our tool (see the Results section)!
Mortality, Replication, and Saturation
The first thing we know about these populations is that they can die after a certain period and replicate at a certain rate. For safety reasons, the bacteria in population \( C \) carry a dapA mutation and are unable to replicate in the absence of diaminopimelic acid, which they won't find in the environment. Therefore, there is no replication term for population \( C \). Population growth is represented by the coefficients \( r_i \).
These populations exist in an environment containing a certain amount of energy and nutritional resources. As a result, the populations cannot increase exponentially to infinity. Self-regulatory phenomena will, therefore, come into play. Hence, our model will be a saturated population model:
Verhulst model and Allee effect model
These phenomena are taken into account in the Verhulst (1838) model, also known as the logistic growth model. This continuous-time model takes the form of a differential equation:
With "r" being the population growth rate and "K" being the carrying capacity of the environment, which is the maximum number of individuals the environment can support considering factors like space and resources. A logistic growth population always tends toward "K," the carrying capacity, regardless of the initial population density (\( P0>0 \)). This implies that as long as a population is not entirely eliminated, it will continue to grow toward the carrying capacity "K." This means, for example, in the case of resistant bacteria, as long as we haven't completely eradicated them, they will persist in reproducing and occupying available space. Another way to model the carrying capacity of the environment is by adding an Allee effect to the logistic model. Incorporating an Allee effect accounts for the fact that below a certain threshold, the population will rapidly decline toward zero. The Allee effect is found in many population dynamics models. This effect can be modeled as follows:
" \( A \)" corresponds to the critical density below which the population growth rate becomes negative. In our model, "\( A \)" is not tied to the population density of resistant bacteria but rather to the ratio of resistant bacteria to the total population. This is because below a certain threshold in this ratio, the action of our tool is considered successful, as the population of resistant bacteria is destined to decline. In a scenario like the human body, this could imply that the immune system has taken over. We model the Allee effect only for populations of resistant bacteria since these populations could otherwise diverge from zero due to the combination of resistances. We incorporate this Allee effect into both the replication and conjugation terms of resistant bacteria growth. The ratio is defined as follows:
Therefore, we will model according to the model with the Allee effect, as it best aligns with the biological situation, taking into account the limiting resources, the ratio of sensitive to resistant bacteria, and in the case of the human body, the metabolic processes that can take over to eliminate a small quantity of resistant bacteria.
Passage between compartments
The populations are represented as compartments, and bacteria in one compartment can transition to another if it is possible and according to certain parameters. Therefore, we are approaching an epidemiological model of the SIRS type. Indeed, our model is based on this type of model because it adheres to its two fundamental principles: compartments and rules. Compartments categorize the population into various possible states concerning the disease, while rules specify the proportion of individuals moving from one compartment to another. The added complexity compared to a typical SIRS model is that our populations vary (due to mortality and "births"). Resistant bacteria can acquire the CRISPR tool plasmid: \( R \)->\( Rc \), \( S \)->\( Sc \). The plasmid used has low stability and can be lost at any time, causing bacteria in the \( Rc \) population to shift to \( R \) and those in the \( Sc \) population to shift to \( S \) (as well as \( C \) to \( S \), a transition we will neglect for now). Once induced transcription occurs (after the introduction of anhydrotetracycline into the environment), resistant bacteria containing the tool plasmid become sensitive (we do not consider the transition \( Rc \)->\( Sc \)->\( S \) here; we assume that after transcription, the plasmid is lost, leading to the transition \( Rc \)->\( S \)). Furthermore, resistant bacteria carry this resistance on a conjugative plasmid. Hence, sensitive bacteria can then become resistant \( S \)->\( R \) and \( Sc \)->\( Rc \). We assume that the reverse transitions do not occur because the plasmid containing the resistance is stable (if it were not the case, we would face the global problem we are trying to combat today: the spread of strong and enduring resistance over time). Bacteria in the \( C(t) \) population can be resistant or sensitive, and we will neglect the difference since they will lose resistance after induction; we consider them as sensitive here. This yields the figure below:
For the transition of resistant bacteria Rc to sensitive, as well as sensitive \( Sc \), and \( C \) to \( S \), which occurs after induction, we require a specific term. This term will depend on when we induce the transcription of our tool plasmid. To ensure the self-contained nature of our system of differential equations, we have introduced a sixth variable: \( A(t) \), for anhydrotetracycline. For all transitions, the terms consist of a constant representing the specific characteristics of that transition, multiplied by the involved variables. This holds true except for the transition after induction of all bacteria possessing the tool into the \( S \) population. In this case, we introduce a function, theta, which applies to the variable \( A \). Theta will yield a zero value until \( A \) is sufficiently large, and then the function will increase as \( A \) becomes sufficiently large according to a Hill equation. We have chosen to use a Hill equation so that the transitions to \( S \) are smoother, controllable through the parameters of this function, and mathematically modeling makes biological sense. Hill equations were developed based on biological processes, making them more relevant in our context than, for example, an \( e^{-x} \) function, which lacks biological significance. The theta function is as follows:
Addition of anhydrotetracycline to the medium
At the beginning of our study, at time \( t=0 \), there is no inducer in the environment. We select a time, \( t_1 \), after the \( C \) bacteria have had the opportunity to transfer their plasmid to a significant, if not all, portion of the bacteria through conjugation. At this time \( t_1 \), we will alter the initial condition of \( A \) to introduce the inducer into the environment. This process follows the following model:
The value of variable \( X_1 \) at \( t_1^+ \) will be updated based on the value at \( t_1^- \) applied to the function \( g \). In our case, the function \( g \) will be a constant function returning a specific quantity of anhydrotetracycline.
Conjugation efficiency and encounter probability
The efficiency of conjugation will depend primarily on two factors: the conjugative power of our plasmid and the probability of encounter between a donor and a recipient. For the time being, we are assuming that conjugation occurs instantaneously. However, it is likely that we will need to introduce a delay into the model later on. Conjugation requires a donor to encounter a recipient. The challenge here is that this takes us down to the bacterial scale, whereas the model operates at the population level. We can incorporate this probability into an individual-centered model. For now, we will retain a term dependent on the number of resistant bacteria and the number of donor bacteria. We can attempt to give more weight to the number of donor bacteria, for example, by introducing a quadratic term.
We have therefore at this stage and after numerous discussions, revisions and corrections of the model, designed our system of differential equations. Before going directly to the simulations, we will check that the model is likely to tend towards the result we want, that is to say: only sensitive ones at the end of the experiment. To do this, we verify by theoretical analysis (and therefore by hand) that there is an equilibrium point corresponding to this situation. An equilibrium point in a system of differential equations is a particular state where the rates of change of all variables in the system are zero. In other words, at the equilibrium point, the variables in the system do not change over time. In mathematical terms, an equilibrium point is reached when all derivatives of the system variables with respect to time are equal to zero. This means that the variables no longer vary, and the system is stable in this state. We therefore carry out an analysis of the equilibirum points :
We therefore obtain a balance point corresponding to what we want :\( R^{*}=0 \), \( R_c^{*}=0 \), \( S^{*}=K \), \( S_c^{*}=0 \), \( C^{*}=0 \). Its existence is guaranteed here. However, it is also necessary to check that this equilibrium point is stable (we will actually reach it). To check the stability of our balance point, manual calculations are impossible (here we have a matrix of 36 elements to calculate, we generally stop at 9 in manual analysis...). We will therefore move on to the simulations. We have no assurance that we will indeed achieve balance in a stable manner, but we already know that this balance exists, which is very reassuring for the future.
Choice of parameters and initial conditions
Our bacterial strains are E.Coli. Under optimal conditions (37°C, nutrient-rich medium, no antibiotics, good oxygen content, etc.), sensitive bacteria divide every 20 minutes or so. Resistant bacteria containing a sizeable plasmid will take longer to divide, so we'll choose a population doubling time of around 30 minutes. The impact of tool plasmid ownership on doubling time will be neglected here, to facilitate the model.
To calculate the growth rate of the populations of susceptible and resistant bacteria from the doubling times, we will use the exponential population growth model:
We can then obtain a link between r and the doubling time td as follows. The values obtained for \( r_R = r_{R_c} \) and \( r_S = r_{S_c} \) are summarized in Table.1 :
The same applies to\( k_A \), the decay parameter for anhydrotetracycline, which has a half-life of around 6 hours.
Our model is now complete, and we have all our parameter values, so we can run our simulations. The python code used is available on our gitlab. The equations are approximated using Euler's method. The following graphs show the evolution of \( C, Rc, S, Sc, R, ARS \) and \( ASR \) as a function of time. At time \( t=2h \), there is induction by anhydrotetracycline.
The results are conclusive regarding the effectiveness of our tool. At time \( t=0 \), we have \( R \) bacteria, \( S \) bacteria, and \( C \) bacteria. The \( R \) and \( S \) bacteria will mostly become \( R_c \) and \( S_c\). Before induction, there are transitions between different populations due to plasmid losses or resistance conjugation. When we induce our tool plasmids, all the bacteria that possessed them become \( S \). The population of sensitive bacteria then becomes significantly dominant, leading to the extinction of resistant bacteria. Thus, we observe a pronounced growth curve for the \( S \) population.
The in silico study of our CRISPR dCas9 plasmid tool concludes that it works!
The modeling program built based on our system of differential equations allows us to test various initial conditions. Therefore, prior to conducting experimental tests, laboratory team members can anticipate the quantities of resistant, sensitive, and tool bacteria they will introduce into the environment. Other parameters that can be adjusted to predict population behavior include the inducer concentration added to the environment and the timing of the induction of tool plasmid transcription.
Our model is therefore a major support for predicting the behaviour of bacterial populations in the laboratory testing phase!
We have developed a version of our visual model suitable for people suffering from color blindness (more informations on the Inclusivity page).
Double-click on the button below for the inclusive version!
In order to visually illustrate the evolution of our bacterial populations in a more accessible manner, we employ a model inspired by individual-based models (IBMs). IBM models were initially used in ecology and are based on an explicit representation of all individuals in the system, in contrast to dynamic models. They leverage individual behavior within a group to determine individual and collective responses to changes in the environment, such as alterations in food conditions or increased disturbances. In our case, we will represent our bacteria as spheres moving in a 2D plane. The size of the plane, initial positions of the bacteria, as well as their speeds and directions of movement will be entirely arbitrary. The objective here is not to obtain actionable results but to present the model of equations in a more understandable way for non-mathematicians. The bacterial populations are represented as follows:
The bacteria are individually represented in a 2D space. We use the same initial conditions as in our differential equations model: 5 dark green balls (\( S\)), 5 red balls (\( R \)), and 10 blue balls (\( C \)). The bacteria will move within the plane. The growth, mortality, conjugation, and plasmid loss parameters are the same as those used in our differential equations model. We summarize the possible interactions and evolutions as follows:
We will allow the bacteria to move and evolve for 2 hours, and at \( t=2 \), we will induce anhydrotetracycline. For simplification, this induction will involve converting all light green, orange, and blue balls directly to dark green. The final frame of the resulting video shows the experiment's outcome: if all the balls are dark green, the experiment was successful. This simplification remains relevant because, in the model, the few remaining resistant bacteria (or Rc bacteria that might have lost the plasmid during induction) would have disappeared due to the Allee effect. The simulations involve an element of randomness and are never the same when launched. However, across all simulations run, the result is always the same: the experiment is successful. To best illustrate our model, we started by conducting a simulation with only dark green (sensitive) balls and red (resistant) balls at the beginning. This simulation serves as an illustration of what happens when our donor bacteria are not introduced into the environment. It can be observed that very quickly, the sensitive bacteria acquire the resistance plasmid, and the environment becomes populated solely by resistant bacteria.
We then run simulations of our complete IBM model. Here's an example of the resulting simulation:
This model does not provide additional results or more in-depth analyses. However, it offers an extremely valuable asset: accessibility for non-mathematicians. Indeed, our team is committed to being multidisciplinary, and each member aims to grasp all aspects of the SuperBugBuster project. This effort to simplify complex concepts is also undertaken to ensure that anyone can understand the work presented on this page. Effective and conclusive work relies on communication between biologists and modelers.
The code for this model, like that of the differential equations model, can be found on our gitlab.
One of the values of iGEM is cooperation between teams. So we've decided to include on our model page a few tips and rules to help future iGEM teams model their projects. We'll explain some of the main types of equations and how they can be used in different situations.
Our model was constructed based on several simplification assumptions. We justified these assumptions, as the model retains biological relevance, but it could be further complexified to better correspond to biological reality. The complexifications we have considered so far include:
The envisioned implementation of SuperBugBuster was to use it as a probiotic for an effect on the intestinal flora. The intestinal flora is composed of a multitude of different bacterial populations. Therefore, an improvement to the model could involve adding different variables \( S_i \) corresponding to the populations of bacteria, each with its own unique characteristics.
https://www.apmep.fr/IMG/pdf/AAA05047.pdf
NetLogo Home Page. https://ccl.northwestern.edu/netlogo/
« Modèles compartimentaux en épidémiologie ». Wikipédia, 6 septembre 2023. Wikipedia, https://fr.wikipedia.org/w/index.php?title=Mod%C3%A8les_compartimentaux_en_%C3%A9pid%C3%A9miologie&oldid=207634569
« Effet Allee ». Wikipédia, 1 septembre 2023. Wikipedia, https://fr.wikipedia.org/w/index.php?title=Effet_Allee&oldid=207470569.
Jaumain, Roland. Utilisation de la conjugaison bactérienne, associée à Crispr/Cas9 pour combattre les bactéries multi-résistantes. 2021.
Écologie / Ecology - Modèle individu centré. https://www.bonobosworld.org/fr/glossaire/modele-individu-centre
