Overview of current models
Deterministic models
The modelling of plasmids started with the work of Stewart and Levin in 1977 studyng the stability on conjugative plasmids using a mathematical model of ordinary differential equations (ODEs). For this reason we will begin to analyze deterministic model involving conjugation.
Deterministic models regarding conjugation
The first model of ODEs for the study of conjugation was developed in the work presented in \citep{nokey_b7e57}. The model became popular and utilized in posterior studies because it is very complete in the processes affecting the phenomena. It consists on two populations: one bearing a plasmid and another plasmid free, each population has a hyperbolic growth function proposed by Monod, the burden of bearing a plasmid is accounted by differences in parameters of the growth function and the quantities of resource required to cell division, the segregational loss is incorporated as a rate and conjugation by a rate and it also accounts for selection by as fitness cost. The model then considers two distinct regimes of resource utilization, one called "equable" simulating a chemostat setup and other "seasonal mode" for modelling batch cultures. The major conclusion of the work is that "There is a broad range of condition where conjugative plasmids can become established and where plasmid-bearing bacteria will maintain high frequencies, even when these factors considerably reduce the fitness of their host".
They continued their work in a simpler model postulated in \citep{Levin_1979}. This model is constructed under the following assumptions: conjugation occurs at a rate that is proportional to the concentration of plasmid free and plasmid bearing cells, there is no segregational rate, a transconjugant cells can act immediately as donor, original donors and transconjugants conjugates at the same rate, there is no plasmid cost. The first assumption can be interpreted as the mass action law and its the parameter ruling the dynamics the system. The work concludes that plasmid transmission is well approximated by the mass action model with constant rate of transmission but only during exponential growth. Experimental data show that during lag phase the rate is higher whereas in stationary phase the rate is lower.
In a later work by \citep{Simonsen_1990} the above model was modified to describe plasmid transfer in batch cultures. To do so, the authors added a resource concentration variable and modified growth functions, as well as the conjugation rate, to a Monod function dependent on the resource concentration . Then performing transformation of the ODEs system they found a formula for estimating the transfer rate from experimental observations. They found that the occurrence of a lag phase has no appreciable influence in the estimation of transfer rate, that it also is unaffected by the initial donor/recipient ratio. And that the transfer rate is constant.
So far we have examined studies regarding plasmid stability, in \citep{Svara_2011}, the authors explore the selective pressures which favors plasmid-carried resistance genes, such as dosage of antibiotics and interval between treatments. To address this issue they postulate a ODEs model under the next assumptions: logistic growth with a growth rate, there is a density dependent cell death rate, selection by antibiotics acts by fitness cost, there is no compensation, conjugation occurs at a constant rate and resistance is carried in the chromosome or a plasmid. Based on analytical and numerical results they concluded that transmission is the most important factor promoting plasmid resistance. and that plasmid carried resistance is favored under more intense but less frequent treatment regimes. And a very interesting remark "plasmid-coded resistance will outcompete chromosomally-coded resistance under long-interval treatment regimes".
There is a work involving the relation of cooperation and mobile plasmids presented in \citep{Nogueira_2009}. The aim of this study is to address the question of how can cooperation persist? When the production of public goods is often costly. The work is based in that increased local relatedness favors cooperation. To address this fact they postulated a model using a standard recursion equation for relatedness in a patch-structured population allowing horizontal gene transfer at a constant rate and gene loss at a constant rate. An initial result postulates that invasion of cheaters in a population of cooperators would be prevented if the social trait were coded in a conjugative plasmid by re-acquisition of the trait. As relatedness increases because horizontal gene transfer they suggest that cooperative traits should be preferably coded in the most mobilizable regions of genomes, which is the case for secreted and outer membrane proteins. The major conclusion of the work is that horizontal transfer promotes cooperation.
Deterministic models regarding non-conjugative plasmids
The first effort for modelling non-conjugative plasmids is presented in \citep{Levin_1980}. They present a model for the population biology for non-conjugative plasmids that can be mobilized by conjugative factors. The model is an extension of the one used in \citep{nokey_b7e57} restricted to the chemostat version. They consider four populations: plasmid free (receivers), donors carrying non-conjugative plasmids, mobilizers carrying the conjugative plasmid and mobilized donors which carry both plasmids; each population growth rates are modeled by a linear functions dependent on resource; a unique conjugation constant rate; segregational loss occurs at a constant rate and mobilized donors lose single factor at that double rate. With such a model they conclude that "it is highly unlikely that even readily mobilizable non-conjugative plasmids would become established and could be maintained in bacterial populations in the absence of some form of direct selection favoring their carriage".
The previous work provided insights for the stability of plasmids but with the influence of a conjugative plasmid. In the next publications we show simpler models for the dynamics of a population carrying a single non-conjugative plasmids.
We will start with \citep{Lenski_1987}, the work consists in postulating a single ODE for the frequency of plasmid-bearing cells that only takes into account the rate of segregation and the selective disadvantage (cost) of bearing the plasmid. Integrating the differential equation they get an expression from which they state that the effect of selection is higher when the initial plasmid frequency is low. Then they parameterize the model using a culture of E. coli bearing the plasmid pACYC184 growing in a chemostat and plating samples on selective media. And found that the best fit is reached when segregation rate is 0 and selection coefficient is 0.304. With this observation rather than state a relevance of the parameters they say that segregation and selection must be considered simultaneously when one is addressing the causes of plasmid instability.
In the same year in the study presented in \citep{Cooper_1987} the authors postulate a model similar to the previous but for estimating the fraction of plasmid free bacteria. The parameters of this model are the segregation rate and the difference of growth rates between the plasmid free and the plasmid bearer bacteria. Then by making assumptions of the relation of this two parameters, they reduce the model to a three case specific models and compare the analytical results to experimental data obtained from bacterial growing in continuous cultures that resemble the assumption, for example for growth rate difference much grater than segregation rate they used E. coli 1B373 with the plasmid pMG169 which has a cost of 3.06x10-2 and a segregation rate of 8.1x10-6 and successfully fitted the data. Except from a misbehaved part of a experimental data set which they propose, could be due to compensatory mutations. This could be consider the major contribution of the work, and one case of the model could be used to fit this part.
In the another work described in \citep{Stephanopoulos_1988} the authors aim to study stability of of recombinant (plasmid bearer) cells growing in a chemostat by analyzing stability of an ODEs model using theory of index of a singular point. The model is constructed under the next assumptions: a culture consists of two types of cells parental and negative, parental cells produce negative cells variants at a constant probability, negative variants cannot revert to give parental cells, parental and negative variants may have different growth rates. Then the model consist in three ODE, one for cell variant and another for substrate concentration with growing functions substrate dependent. The stability analysis is then performed by finding the contour curves when each of the population respective ODE is set to be equal to 0. By doing it they found that a necessary condition for stability of the coexistence steady state is that the parental cells grow faster than the reverant cells in, at least some region of the limiting substrate concentration. Which in general is not expected unless plasmid compensatory mutation appear.
In a study involving plasmids with different copy numbers presented in \citep{Ayala_Sanmartin_1989}. They propose a simple ODE model for the stability of plasmids assuming random segregation. And experimentally determined the stability and plasmid copy number (pT7: 11, pBR322: 61, pIO1: 37, pBR327: 74, pBR329: 20) as well as the degree of multimerization (pT7: monomer, pBR322: dimer, pIO1: monomer, pBR327: dimer, pBR329: monomer). The model satisfactorily fit the data for mid copy number plasmids but failed to predict the stability of high copy number plasmids. They only venture to conclude that multimerization is related to copy number possibly due a saturation-like process.
The greatest motivation of this models to study the stability of plasmids under chemostats for it is a major concern in the manufacture of a product in biotechnology. Such is the case in the next research presented in \citep{Yuan_2010} a completely theoretical work focused in studying the stability on a competition model between plasmid bearing and plasmid free bacteria growing in a pulsed chemostat with a plasmid produced an inhibitor of the more competitive plasmid free cells. The model consist on four ODEs: one for each population, one for the concentration of inhibitor and one for the concentration of substrate. They then make assumptions over the relationship between the maximal growth rate and the dilution rate for each population to find the conditions necessary for the system to oscillate, that is, the conditions for the stable coexistence of the two populations. And conclude that in a chemostat this conditions can be achieved by controlling the initial concentration of substrate, the dilution rate and the period of the impulsive effects.
The above model postulates a way for the coexistence of plasmid bearer and plasmid free bacteria using inhibitor, namely an antibiotic. In the project presented in \citep{Yurtsev_2013} the authors explore the population dynamics of a co-culture of bacteria, one carrying a multicopy plasmid encoding for a b-lactamase enzyme and the other plasmid free. The authors then measure the changes in the fraction of resistant and sensitive bacteria with respect of a initial fraction over a 1-day growth on LB media containing a 50-fold sensitive bacteria ampicillin MIC. The data showed that the population reaches and equilibrium (~0.25 of resistant bacteria). Then used this data to generate a 'difference equation map' and postulated a model that fitted the data. The model takes into account growth rates, sensitive cells death rates, antibiotic concentration and degradation, initial fraction, final fraction, cell densities, and a sub-population specific parameter for the duration of the lag phase. By analytically solving model they get a formula for estimating the equilibrium-resistant fraction that depends only on the initial cell density and antibiotic concentration. They found, as expected and proved experimentally, that the equilibrium fraction increases as the antibiotic concentration increases. They argue that compensatory mutation should not significantly change the equilibrium fraction but only the time to settle. And that using a high concentration of antibiotic that gets to kill resistant cells will only increase cooperativity.
The previous models care about stability of plasmids but none of them pay special attention to the segregation dynamics which can affect the rate of generation of segregants and thus the stability. The work exposed in \cite{M_nch_2019} engage plasmid stability by focusing on segregation. The approach used is to identify "evolutionary stable strategies", that is the situation that maximizes the average fitness of the population. The model consists in two levels: plasmid reproduction within bacteria and bacterial division. Both processes are modeled by a logistic function with respective growth rates parameters. When a bacteria divides, the segregation process is describes by a binomial distribution B(N,p) with N being the plasmid copy number and p (probability of success) being the symmetry of segregation p=1/2 for equal (random) segregation. By modifying the latter parameter the model account for unequal segregation. The major finding of this project is that for low copy plasmid bacteria present equal segregation while high copy plasmids one of the daughter cell distinctively and systematically receive more plasmids.
In the next theoretical work \citep{Kentzoglanakis_2013}. The authors apply the perspective of policing and obedience plasmid copy number control mechanism to plasmids R1 and ColE1 by calling the policing agent antisense RNA and obedience its binding affinity. They postulate a two level model in which the population growth is given by an ODE that incorporates, plasmid cost, plasmid copy number and plasmid fitness contribution to the cell, and plasmid replication rate is fashioned by an equation that incorporates a basal replication rate, inhibitor concentration and the affinity of the inhibitor. With this model they first get biologically characteristic copy numbers and find a bounded range of a positive correlation between the plasmid replication rate and the affinity of the inhibitor for the population to exist and maintain the plasmid. That is, for the plasmid to be stable.
Deterministic models regarding compensation
Compensatory mutations are another often ignored in modelling plasmid stability but they could be key elements in explaining unexpected behaviors in plasmids dynamics. As is the case presented in \citep{Millan_2014}, where they studied the stability of a costly non-conjugative plasmid pNUK73 in a population of Pseudomonas aeruginosa under non selective media. From the experimental stability assay the authors found that after a period of exponential lose, the plasmid stabilizes. Contradicting their theoretical predictions using an ODE model based on a random segregation rate, plasmid cost, growth functions depending on substrate concentration. At the end of experiment they performed competition experiments and find that the plasmid cost had substantially decreased. So they modified the model to integrate compensatory mutations both in plasmid bearing and plasmid free populations and explained that the observed stabilization of the plasmid is the joint effect of the decrease of the plasmid bearing population and the slow increase of a compensated plasmid bearing population. This fact was latter confirmed by sequencing, and also found that compensation were due only to chromosomal mutations.
Under the light that compensatory mutations could occur either in the chromosome or the plasmid, the authors in \citep{Zwanzig_2019} explore the different plasmid dynamics that genome localization of compensatory mutation can entail. For chromosomal mutations are inherited vertically and plasmid mutation are spread both vertically and horizontally. To asset this issue they developed an ODE model under the next assumptions: the system is well-mixed, bacteria grow with a maximal growth rate, plasmids entails a cost, compensatory mutations ameliorate that cost and occur a constant rate, bacteria are lost by dilution rate and antibiotic killing rate, there is constant segregation rate, and plasmids are transferred horizontally at a constant rate. They performed simulations to show that compensatory adaptation allow plasmid mediated antibiotic resistance to persist even in the absence of antibiotics and that plasmid mutations facilitate plasmid persistence even when the amelioration effect is is lesser than the provided by chromosomal mutations.
Deterministic models regarding PSK
The next theoretical approach the address the dynamics of conjugative plasmids carrying antibiotic resistant genes incorporating many threats as conjugation and a PSK system in a chemostat scheme with antibiotics cycling \citep{Willms_2006} and focused in the role of recipient cells on plasmid stability. The model is based in that postulated in \citep{nokey_b7e57} adding the effect of a bacteriostatic antibiotic, for which plasmid free cells are susceptible and plasmid bearer cells are totally resistant and a PSK system through the effect of growth cost when a missegregation event occurs. By performing simulations they found that the resource concentration and the PSK efficiency allows substantial gains in plasmid persistence and the greatest effect is achieved when resource is moderate to high and high PSK efficiency. They also state that when cycling antibiotics, recipient cells is the more important factor affecting plasmid frequency.
Further analysis in PSK systems is boarded in \citep{Rankin_2012} where by a mathematical model, they examine the conditions driving the evolution and evolutionary stability of plasmid toxin-antitoxin systems. The model is a susceptible-infected type with three populations: one for plasmid free cells, other for a plasmid bearer population that does not posses PSK system and another population carrying a plasmid with the PSK system, as the previous model, they introduce a term for the effect of the PSK. They also do not allow coinfection, and allow integration of the PSK system to the chromosome, conjugation and segregational loss. With this model they explore the emergence of PSK systems on plasmids from two views, toxin-first and antitoxin-first. And conclude that either way the trait must offer and advantage in at least one environment. They also used a rock-paper-scissors approach to help explain the presence of toxin-antitoxin genes in chromosome.
Stochastic Models
This models uses stochastic differential equations (SDEs) are differential equations that incorporates a stochastic variable or a term that typically incorporates noise to the system.
The use of SDEs in the study of plasmids stability is based in the following project. In \citep{De_Gelder_2004} they monitored the loss antibiotic resistance in a population E. coli K12 carrying the multiresistance IncP1-B plasmid pB10 during a ~500 generation batch experiment without selection. The plasmid is a 64.5 kb conjugative plasmid, with from data from the same experiment the plasmid presents high maintenance. Nonetheless, they observed that the loss of tetracycline (Tc) resistance increased slowly from 0.1 to 7%. To understand this process they postulated an ODE model for the fraction of mutant plasmid bearer population with the following assumptions: mutations are non-reversible, so cells can lose Tc resistance but cannot gain it; mutants increase by growth of previous mutants and from mutations on wild-type cells. With this simple model they successfully predict the observed data and helped to explain that the dynamic observed in the population is due to a high mutation rate and low but significant selection coefficient.
Later in a work from the same group \cite{Ponciano_2007}this model was used to build an stochastic model and compare it to the deterministic model in fitting a time-series data on plasmid persistence on seven bacterial strains. In the deterministic case, they incorporate conjugation on \citep{Levin_1979} but relaxing the mass action assumption and multiplying conjugation rate for a Michaelis-Menten dynamic instead of population size. For the formulation of the stochastic model they considered that there could be periods of heavy plasmid loss followed by stable segregants frequency. Incorporation this to the previous model they postulated that the solution to the system becomes a Markov process with respective probability density function (pdf). With the models they successfully fit three out of seven data sets. Concluding first that different strains presents different plasmid loss dynamics. That the plasmid cost changes form one day to the other, so modelling plasmid cost as random effect substantially improves their ability to adequately explain much of the data.
Another model involving modelling conjugation with stochastic differential equations is presented in \citep{Philipsen_2010}. In this work they analyze data from an experiment for conjugation between two Enterococcus faecium species growing in exhaustible media. The contribution of this work is taking into account that conjugation rate is substrate dependent. The cost of plasmids is implemented as a constant rate that is multiplied to a substrate dependent Monod growth function. The stochasticity is implemented by an additive noise function that is based on noise the observation of the experimental data. At the end they do not make any biological conclusion other than the importance of taking into account that conjugation rate is substrate dependent which best describe observations of conjugation during lag and stationary phases.
In a later work , with aim to investigate the fixation probability of a gene that can be horizontally transferred. The authors in \citep{Tazzyman_2013} postulates two different modelling approaches a branching process and a diffusion approximation expecting to avoid results that are artifact of the model rather than the process under study. The initial model assumes a fixed population size of bacteria, bacteria life cycles has different period of growth and division and mutant population with a low initial frequency. So work consists in finding the fixation probability of a mutant gene that can be horizontally transferred which was never done using diffusion approximation. They found a a previously demonstrated fact, that there is trade-off between horizontal and vertical transmission. That the fixation probability of a deleterious allele can be non-zero, in other words they show a way in which antibiotic resistance genes could become fixed in a population even in the absence of antibiotics.
Wright-Fisher Model
The Wright-fisher model for genetic drift that provides powerful insights into population genetic dynamics. The model is not popular because is build with very specific assumptions that are not often satisfied. This assumptions are a fixed population size, there is no selection, no mutation, no migration and non-overlapping generation times and diploid populations.
In the study presented in \citep{Ilhan_2018} they apply a standard haploid version of the Wright-Fisher model to simulate the evolution of a population that is subject to random genetic drift \citep{guillespie2010}. To incorporate plasmids evolution they followed the approach used in \citep{Peng_2005} a study of mitochondrial evolution. With this framework they study plasmid segregational drift, a process well described in mitochondria but never studied before in plasmids. Using both modelling and experimental approaches they compared the relation between plasmid copy number and mutations (or allele variants) drift. Comparing populations with different plasmid copy numbers and chromosomal alleles (population sizes adjusted). They conclude that plasmid occurring mutations are easily lost by segregational drift and that the allele frequency of mutation residing on plasmids is increasing slower in comparison to chromosomal alleles.
Individual Based Models
Individual-based models (IbMs) are based in the view of the uniqueness of individuals in a system. In this models each individual has its own state of a set variables which depending on the phenomena under study can undergo a very complex computationally. This complexity relies in that simulation are heterogeneous in processes and interactions, and therefore increase the chances of programming errors and testing/validation time \citep{Gregory_2006}.
The ecological modelling with IbMs started with HERBY \citep{Devine_1997} a discrete grid based environment of plants with herbivores populating the grid, the algorithm includes learning and moving costs and can be applied to sexual and not sexual populations. The next step was the program COSMIC \citep{Gregory_2004} which aim was to simulate bacterial evolution in a multi-scale perspective in the sense of incorporating genes (and mutations), gene products and thus, individuality which in time was affected by environment. At the same time RUBAM \citep{Vlachos_2004} was developed with the aim to simulate adaptive behavior using a grid based environments with nutrient gradients and introduced antibiotics. Than integrating both COMSMIC and RUBAM strategies came COSMIC-Rules \citep{Gregory_2006} this program works at three organizational levels : genes, cells and environment, the last one includes other cells and therefore could be used to study bacteriophages or plasmids.
COSMIC-Rules was used in a completely theoretical work to study the transmission of plasmids in a bacterial population \citep{Gregory_2008}. Their experiments consisted on a comparison between invasions with conjugative plasmids bearer populations: invading with two incompatible plasmids vs invading with two compatible plasmids, each plasmid carrying a different antibiotic resistance gene and then add both antibiotics to the environment and analyze plasmid dynamics. In this work conjugation is based on proximity and metabolic state, the latter meaning that transconjugant cells have a waiting time before they can act as plasmid donors and so the donor cell before donating to a second recipient. Plasmid cost was taken into account as 1% growth rate and were not specific for double transconjugants. Selection with antibiotics was boolean and assumed 100% vertical transmission. it also incorporates a parameter for surface exclusion so the same cell can not act as receiver of the same plasmid type to model plasmid incompatibility.
In this spatially structured study they found that the optimal conjugation rate was 1x10-3, which is rather high. From incompatible plasmids invasion experiment they concluded that this property is a limiting factor for the spread of plasmids. From the compatible plasmids experiment they concluded that both are able to be transferred and maintained in the entire population. The major contribution is to state a prove of concept of this approach to study gene transfer by conjugation.
The above model presents a general overview of plasmid dynamic behavior with major plasmid properties, a more property-specific model is proposed in \citep{Merkey_2011}, with the motivation of understanding poor plasmid spread in deeper biofilm layers the authors postulate an IbM consisting in three agents: cells, extracellular polymeric substances (EPS) and plasmids. Considering plasmids as another agent they can simulate different plasmid types. They also take into account plasmid burdens and metabolic reactions, a segregation probability, conjugation, conjugation lag times and conjugation specific parameters such as pilus reach distance and pilus scan speed. They also introduce the concept of growth dependence of HGT throughout the pilus scan speed parameter modeled as piece-wise linear function dependent a growth tone parameter with specific cut-offs. With this model they explore the influence of the growth tone during invasion experiments and were "able to reproduce observed dynamics of plasmids in biofilms, including cases of complete invasion and invasion limited to the biofilm surface". They also found that timing parameters, distance between neighbors as well as plasmid burden have a stronger influence in the invasion success than the growth rate of the receivers or the segregational loss rate. The relevance of this works relies on the insights of explaining the limited spread of conjugative plasmids in some biofilm communities by the growth dependence of the HGT process.
Latter in a work by the group of Michael Brockhurst \citep{Harrison_2016} they postulate an enigmatic IbM with the aim to simulate the dynamics of the pQBR103 mega plasmid that carries a mercury resistance cassette within a transposon and a fitness cost of 25%. With their model they explore a mutation rate parameter that compensates the fitness cost to 99% and transposition rates and conclude that the stability of plasmids is a race between incorporation of the accessory genes into the chromosome and the appearance of compensatory mutation. And that in the absence of selection mutations must have large fitness effects and occur at high rate.
In the search for resolving the plasmid paradox, a great contribution was made in the work made by \citep{Werisch_2017}. Using an IbM they propose a mechanism to maintain non-transmissible plasmids. The model contemplates two types of plasmids non-transmissible and transmissible, the last ones with two states: able to transfer or not. Segregational loss and incompatibility only occurs when the two types of plasmids are present in the same cell, and is modeled as per-type segregation probability, accounting that two incompatible plasmid types are unable to persist in the same cell due to cross reactions between replication control mechanisms. They resume their finding in the following : "non-transmissible plasmids can be maintained along with co-occurring, incompatible conjugative plasmids, although non-transmissible plasmids provide no advantage to the host and will be lost in the absence of this plasmid composition".
The last model to analyze is namely an individual-based model in the sense that it is constructed based on individual assumptions but consist in a walkthrough all the modelling approaches revised. In the work presented in \citep{Billiard_2016}, the authors consider two clonal populations of haploid individuals each carrying a different 'trait', which can be interpreted a conjugative plasmid. They also consider natural death of individuals, competition between trait subpopulations. The model is constructed under general horizontal transfer without the possibility of coexistence of both traits, and thus, there is conversion of individuals occurring at different rates. By considering that horizontal transfer is stochastic they formulate an equation for an infinitesimal generator of a Markov process in which each populations are scaled by a specific parameter K, when K tends to infinity, the sequence of the stochastic process converges to the solution of a system of ODEs similar to a Lotka-Volterra. By analyzing the phase diagrams of this ODEs system they find they could be stable coexistence of the traits. Then to see magnitude of the fluctuations around the deterministic dynamics they construct an Ornstein-Uhlenbeck process to which solution is stochastic differential equation with Brownian motions related to the birth death processes of each population and to the horizontal transfer, then using some assumptions about later processes they construct a simpler stochastic differential equation that latter rewrite it as the Wright-Fisher diffusion approximation.
Then they use this last model to analyze the invasion (and fixation) of plasmids. The conclusions of this work can be summarized in that fixation becomes case specific and dependent on the transfer rates and costs, for example, in individuals with small birth rate, fixation is more sensitive to transfer rate than to selective value. At the end the take home message is that "horizontal transfer interacts with ecology (competition) in non-trivial ways. Competition influences individual demographics, and this in in turn affects population size, which feeds back on the dynamics of transfer".