Understanding MRSA clonal competition within a UK hospital; the possible importance of density dependence

Background: Methicillin resistant Staphylococcus aureus (MRSA) bacteria cause serious, often healthcareassociated infections and are frequently highly resistant to diverse antibiotics. Multiple MRSA clonal complexes (CCs) have evolved independently and countries have different prevalent CCs. It is unclear when and why the dominant CC in a region may switch. Methods: We developed a mathematical deterministic model of MRSA CC competing for limited resource. The model distinguishes ‘standard MRSA’ and multidrug resistant sub-populations within each CC, allowing for resistance loss and transfer between same CC bacteria. We first analysed how dynamics of this system depend on growth-rate and resistance-potential differences between CCs, and on their resistance gene accumulation. We then fit the model to capture the longitudinal CC dynamics observed at a single UK hospital, which exemplified the UK-wide switch from mainly CC30 to mainly CC22. Results: We find that within a CC, gain and loss of resistance can allow for co-existence of sensitive and resistant sub-populations. Due to more efficient transfer of resistance at higher CC density, more drug resistance can accumulate in the population of a more prevalent CC. We show how this process of density dependent competition, together with prevalence disruption, could explain the relatively sudden switch from mainly CC30 to mainly CC22 in the UK hospital setting. Alternatively, the observed hospital dynamics could be reproduced by assuming that multidrug resistant CC22 evolved only around 2004. Conclusions: We showed how higher prevalence may advantage a CC by allowing it to acquire antimicrobial resistances more easily. Due to this density dependence in competition, dominance in an area can depend on historic contingencies; the MRSA CC that happened to be first could stay dominant because of its high prevalence advantage. This then could help explain the stability, despite frequent stochastic introductions across borders, of geographic differences in MRSA CC.


Introduction
Staphylococcus aureus are commensal organisms and opportunistic pathogens that can be clustered into related individual lineages or clonal complexes (CC) (Jensen and Lyon, 2009;Deurenberg and Stobberingh, 2008). Some isolates have acquired SCCmec, a genetic element that confers methicillin resistance (MRSA). MRSA infections now cause significant morbidity and mortality around the globe (Grundmann et al., 2006;van Hal et al., 2012;de Kraker et al., 2011;Cassini et al., 2019). MRSA clones are resistant to virtually all beta-lactamases, favoured antimicrobials for treatment and prophylaxis. On top of this, MRSA isolates often carry additional resistances, many of which are encoded on mobile genetic elements (MGEs) (e.g. to aminoglycosides). Resistance to all antibiotic classes has been described in S. aureus, but these additional resistances are variably distributed across MRSA CC.
Additional resistance is predicted to confer a selective advantage in the hospital setting, where commensal organisms will be exposed to antibiotics, even if carriage of MGEs also incurs some small cost. The clonal complexes (CCs) mainly evolve their resistance independently, due to restriction modification systems between CCs, which prevent the entry of "foreign" DNA (Feil et al., 2003;Harris et al., 2010). Movement of resistances on MGEs within a CC is likely to be frequent: it has been seen in individual patients and in animal models (Moore and Lindsay, 2001;McCarthy et al., 2012;Goerke et al., 2004;Stanczak-Mrozek et al., 2015). Interestingly, although the CCs of methicillin susceptible S. aureus are widely distributed spatially, a limited number of successful CC clones dominate the MRSA populations of different countries Cockfield et al., 2007). The reasons for this country level segregation in MRSA clones are unknown (Johnson, 2011), although inter-hospital patient transfer networks, linking MRSA populations within countries, likely play a role (Donker et al., 2010(Donker et al., , 2017. Within countries, MRSA clonal dynamics are relatively stable over years to decades, although shifts in dominant clones have been seen, for example in Singapore and Portugal (Hsu et al., 2015;Aires-de-Sousa et al., 2008;Bal et al., 2016). Also in the UK, there has been a switch from CC30 SCCmecII isolates in the 1990s (EMRSA-16), to CC22 SCCmecIV in the 2000s (EMRSA-15) causing the majority of MRSA infections (Wyllie et al., 2011). This was uniquely captured in detail by a case study from a single UK hospital that collected incidence of infection data, CC type, antibiograms and antibiotic usage over the period of change. It also showed the appearance and loss of another MRSA sequence type (subset of a clone) ST239 between the dominant periods (Knight et al., 2012). The ST239 was reported concurrently in at least one other London hospital , although its role in the competition dynamics of MRSA populations remains unknown. In 2020, CC22 remains dominant in the UK (Donker et al., 2017;Coll et al., 2017).
With the current work we wish to address, in general, possible reasons for the spatial and temporal structuring of CC dominance, and in particular what drove the CC switch within the UK hospital described above, also allowing for the relative stability observed since. For this purpose, we developed a new mathematical model to capture the fundamental processes of competition between CCs. Although previous models of antibiotic resistant bacteria such as MRSA have included competing strains (e.g. D' Agata et al., 2009;Kardaś-Słoma et al., 2011)), our model is distinct in that it also includes the processes of loss and gain by horizontal gene transfer of MGEs.
From classical ecological models we know that in general when competition occurs for a single resource, in a closed system, only the competitor who most efficiently uses this resource will in the end survive (Tilman, 1982;Miller et al., 2005). We therefore expect the CC with the optimal balance of fitness factors, such as growth-rate and antibiotic resistance, to out-compete all others. For this reason, modelling studies have often explicitly addressed the puzzle of why many different strains of bacteria nonetheless co-exist (Krieger and Hill, 2018;Davies et al., 2019). Difference in usage of antibiotics may help explain the difference in dominant MRSA CC between countries, favouring particular CCs as the top-competitors locally. But what then caused the UK switch in CC dominance when antibiotics usage did not change?
Intriguingly, we find that in our model the inclusion of within CC dynamics allows for multiple stable states of the system; competition between CCs may be density dependent. This factor allowed us to reproduce the dynamics of relative CC prevalence values observed within the UK hospital.

Methods
Model. We describe the dynamics of competing MRSA clonal complexes (CCs) in the hospital as a system of differential equations (Box 1). These equations describe the change in densities of each CC as dependent on CC characteristics and on the availability of resource. The model is kept simple in that the resource for which bacteria compete represents susceptible patients as well as all hospital and equipment surface area (fomites) upon which MRSA can survive, i.e. it ignores spatial barriers. Within each CC, we distinguish between bacteria carrying standard MRSA resistance only, and those with elements conferring additional resistance. These extra resistances can be lost, gained by gene acquisition, and transferred when isolates from the same CC meet.

Data
With the above model, we aim to understand the dynamic changes as observed within St George's Healthcare NHS Trust hospital from 1999 to 2009, as published previously (Knight et al., 2012). This exceptional study documented the relative abundance of isolates from different CCs over time among infected patients. Consistent with this study, we here consider ST239 as a CC (in origin this is a variant of CC8 (Diep and Otto, 2008) (no other CC8 was observed)). We include only CC22, CC30 and ST239, since our main interest is in what caused CC22 to take over the dominant prevalence position from CC30, and all other CCs were found

Box 1
Formal description of the model.
Here j specifies the CC (j = 1, 2, …, n, with n the total number of competing CCs). m j is the density of the standard MRSA resistant subpopulation, while r j is the density of the multidrug resistant sub-population of clonal complex j. z represents the density of the resource for which the sub-populations compete. i mj and i rj are the exogenous rates of inflow of the two strains respectively. We additionally define total CC inflow i j = i mj + i rj . b j is the clonal complex growth-rate, and c j a proportional cost to multidrug resistance resulting in lowered growth rate of the more resistant strain (c j < 1). d is the natural death or removal-rate of bacteria, which includes removal from the hospital via patient discharge and cleaning of fomites, a is the additional antibiotics induced death-rate, and k j the proportional decrease in the antibiotics-induced death rate due to resistance (k j < 1). g represents the rate of resistance transfer when the standard MRSA and the multidrug resistant strains meet, while l is the loss-rate of resistance. s is the rate at which the standard MRSA strain mutates to gain multidrug resistance. only in small numbers (added together, other CCs made up ~8% of prevalence from 2006 to 2009, and ~0% before).

Parameter values
All parameters are tailored to the above hospital setting, and are shown in Table 1. Our model simplifies by having just two strain versions per CC, standard MRSA (resistant to methicillin and other penicillins, fluoroquinolone and erythromycin) and multi-drug resistant (e.g. additionally to aminoglycosides, tetracycline, fusidic acid, chloramphenicol, mupirocin, trimethoprim and/or co-trimoxazole). In reality the number of different drug resistances carried by the bacteria is found to vary within the clonal complexes, with CC22 carrying fewest and ST239 most resistances on average (Knight et al., 2012). Our model captures this average difference between CCs by the k j parameters, which set a proportional decrease in antibiotics-induced death rate due to resistance. Based on the level of resistance in the data, we restrict k 239 > k 30 > k 22 . However, since it is poorly known how much of bacterial death is avoided by resistance genes, we fit k j to the prevalence data of the CCs (see below under model analysis).
The growth-rate parameters b j are informed by the in vitro doubling time experiments which showed CC22 isolates growing faster than CC30 isolates, and both of these growing significantly faster than ST239 (Knight et al., 2012). The same study of the CCs in this hospital showed that increased resistance within a complex was associated with minimal cost to the growth-rate, with no detected effect for CC22 at all. As the above mentioned experiments were not accurate to detect small fitness differences however, and since carrying extra DNA should theoretically incur some cost, we have included a slight decrease in growth-rate of c j =2%, 4% and 6% respectively for the more resistant strain versions of CC22, CC30 and ST239 (corresponding to the average resistance levels of these CC).
For the rates of resistance loss, gain by transfer, as well as gain by genetic mutation (l, g and s), we lack estimates, and therefore explored different values. For reasons of parsimony, we presumed no difference among the CCs in these factors. Since we expect genetic mutations conferring resistance to be rare in the given timeframe of 10 years, for our main analyses the mutation rate between levels was set to zero (s = 0).
The mean length of patient stay decreased from about 6 to about 5 days due to policy change in 2005 (Knight et al., 2012). The daily rate of bacterial removal due to patient discharge or death then was 0.17 before and 0.2 afterwards. Although we assume this to be the main cause for removal, we also allow for other causes of bacterial clearance; we let the total per day removal rate, d, change from 0.28 before to 0.32 from 2005 onwards. Within hospital antibiotic usage is high (Versporten et al., 2018), and although it is not well known in how much MRSA removal this results, this factor should be of significant impact within a hospital; we set the antibiotic induced death rate for sensitive MRSA approximately equal to the background removal rate, at a = 0.3 per day. Overall antibiotic use and other infection control practices were stable over this period, except for some increase in mupirocin decolonisation after 2006 (Knight et al., 2012). Since mupirocin appears to have had no impact on competition between CCs (i.e. it occurred after the take-over in dominance by CC22 which we aim to explain) we have left this factor out of the current analysis.
In 2009 (the end of the period of interest), 2.9% of patients admitted to St. George's NHS Hospital were found to be colonised with MRSA (Krebes et al., 2011). Information on colonisation or infection with MRSA at arrival of patients was not available at CC level. One explanation for take-over by CC22 as dominant complex in this hospital would be a greater inflow of CC22 during later years. However, we avoid this trivial explanation by setting the inflow of CC22 and CC30 to be equal and constant over the considered time-period of 10 years, each at ~0.75% of incoming patients (i 30 , i 22 in Table 1). An outbreak of ST239 was documented in a nearby hospital in 2003. We assumed ~2% of incoming patients to carry ST239 during these years, lowering to ~0.4% of incoming patients afterwards (i 239 in Table 1), which enabled reproduction of the ST239 dynamics observed.
The final factor to be quantified is the fraction of high resistance among incoming bacteria (i rj /i j ). Many of the newly admitted MRSA positive patients will be returning patients who were colonised during a previous hospitalisation. We cannot explicitly model this re-admission process due to lack of data. However, we do expect that additional resistance would be subject to a greater rate of loss outside of the hospital, where it is not advantageous. We therefore calculate the fraction of highly resistant MRSA at inflow i rj /i j as hr j /(m j + r j ), where h scales how strongly the current hospital level of resistance determines the level within the inflow. In our baseline scenario we set h = 0.5, but we additionally explored other values for h in the full range of 0 (total loss) to 1 (no loss). Since the source of ST239 was likely mainly from elsewhere rather than from re-admission patients, and since ST239 resistance levels were found to be very high in 2003, the fraction i r239 /i 239 was set to 90% during its years of high inflow-rate.
Fit to prevalence data * We assume a constant number of patients within the hospital, so that the total patient inflow rate equals the outflow rate of ~0.2 per day (as the mean length of stay is about five days). Then infected inflow i equals this total inflow rate of 0.2 multiplied by the proportion of individuals infected at hospital entrance, as stated per CC in the Methods section.

Model analysis
All analyses are performed in R (Team R., 2013). We examine dynamics over time using the function ode() from package deSolve (Soetaert et al., 2010), and equilibrium states, dm j /dt = dr j /dt = 0, using function searchZeros from package nleqslv (Dennis and Schnabel, 1996). We also used the function uniroot.all from package rootSolve (Garrett, 2010) to aid finding all equilibria in case of multiple stable states.
Before attempting to fit the model to the hospital data, we examine the properties of the model by considering simplified settings. First, we consider a single CC and no MRSA inflow into the hospital, i.e. i j = 0. Secondly, we consider two competing CCs in a closed setting. To enable examination of the impact of the starting densities of both CCs, we assume their initial resistance levels to be at equilibrium with their own starting density x, i.e. we solve (r j + dr j /dt)/(r j + dr j /dt + m j + dm j / dt) = r j /(m j + r j ) with m j + r j = x. Since gain and loss of resistance are relatively fast processes in our model, this semi-steady state like assumption of the more resistant population fraction is an acceptable approximation. In a third step, we add external inflow of the CCs (i j >0).

Model fitting
Understanding the basic model properties, we turn to representing the dynamic changes as observed within the UK hospital from 1999 to 2009. To enable comparison of various modelled scenarios with the data, which documented the relative abundance of CCs, we recalculate to relative abundances of modelled CC22, CC30 and ST239, as CC j = (m j + r j )/(m 22 + r 22 + m 30 + r 30 + m 239 + r 239 ).
How relative MRSA CC prevalence values changed before 1999 (t=0) is not known. For reasons of parsimony, for our initial conditions, we set the model to be in equilibrium, i.e. for all CCs at t=0, dm j /dt = dr j /dt = 0. Where multiple possible stable equilibria were obtained, we used each of these separately in the subsequent step, which was to solve the equations for prevalence over time.
We consider two main hypotheses for what caused the switch in CC dominance at this hospital. First, with our primary model version, we consider whether the two known disruptions of the system together, namely temporary high ST239 inflow, and a permanent decrease in the average length of stay, could have allowed for CC22 to take-over from CC30. Furthermore, we aimed to differentiate importance of these disruptions. For this purpose, we used counter-factual scenarios; leaving out either ST239 inflow, or the change in the length of stay, we re-ran our main fitted scenario.
Secondly, we considered the possibility that an additional event occurred, whereby CC22 went through an evolutionary change. In this secondary model version, we assume that CC22 gained additional transferrable resistances on MGEs (r-strain inception) only from 2004, that is we set r 22 (t = 0) = 0 and i r22 = 0 up to 2004, and afterwards i r22 = min(0.01,hr 22 /(m 22 + r 22 )), i.e. at least 1% of CC22 inflow is of r 22 type from 2004.
Besides examining these two main model versions, for each version we also considered alternative settings for several parameters, as indicated above, resulting in eighteen alternative scenarios in total ( Table 2). For each considered scenario a parameter sweep was performed with step sizes of 1% for each of the resistance level parameters k j . We define the optimal fit per scenario as the one minimising the sum of squared differences between the observed relative CC prevalence data and model outcome.
Dominance of CC22 in the UK MRSA population has continued since 2009 (Donker et al., 2017;Coll et al., 2017), the end of the detailed data collection at St George's Healthcare NHS Trust hospital. As a final extension, we use our model framework to explore what characteristics would enable a newly introduced CC to replace CC22 in this UK hospital.

Coexistence of strains within one CC
When a single MRSA CC is first introduced to a hospital the bacteria can multiply, but as free hospital surface space and un-colonised Each scenario was fit to the data by finding the k parameters, denoting the proportional decrease in the antibiotics induced removal-rate for resistant bacteria, for which the sum of squared differences (SSD) between model outcome and data was minimal. For examples of the fitting space, see Appendix Fig. 4. The impact of parameters on model fit was explored by changing one at a time; all other parameters were set at their base values (see Table 1; at baseline h = 0.5, g = 1, l = 0.03 and s = 0). * In the secondary model version with CC22 r-type introduced only from 2004, we did not run the model with alternative values for the mutation rate s, since if s ∕ = 0, a higher resistance carrying element is obtained by CC22 from 1999, and we in effect regain our primary model setting.
patients become scarce, growth will balance with the removal of bacteria, by death (including cleaning) and by patient discharge (Fig. 1A). That is, equilibrium prevalence will be reached. In this equilibrium, isolates with standard MRSA resistance and with multidrug resistance level co-exist. The equilibrium proportion of the standard m-strain and of the multidrug resistant r-strain depends on the relative growth and survival of the distinct strains, but also on the balance between the loss, the transfer and the de novo gain rate of resistance. For example, with all model parameters at baseline value, the more resistant r-strain makes up 93% of the CC population at equilibrium, but if we increase the loss-rate l, the percentile of r-strain bacteria decreases (Fig. 1B). The level of resistance in the CC population in turn affects the overall prevalence the CC can reach; when antibiotics usage a is high enough, the prevented bacterial death outweighs the cost of carrying the element, rendering carriage of the additional resistance a net advantage. In this case a lower population fraction carrying the extra resistance will lower the total CC prevalence (Fig. 1B, black line).
If we assume no de novo resistance acquisition events to take place, i. e. s=0 (at least during a delimited time-period), there is also the possibility of a steady state with only base resistance, i.e. the more highly resistant r-strain might simply never be introduced into the hospital (dotted lines Fig. 1B).

Competition between CCs
Next, we consider what occurs when two CC are present in a hospital. Generally, when competition occurs for a single resource, in a closed system, it is expected that in the end only the competitor who most efficiently uses this resource will survive. Within our model, efficiency is determined by a CCs growth-rate and by antibiotic resistance level. In a simplified version of our model, without resistance gain by horizontal transfer, i.e. g=0, or conversely without loss of resistance, l=0, we indeed see such simple competitive exclusion; if complexes are equal except that one grows slightly faster, i.e. b 1 >b 2 , that complex will eventually replace all others (see Appendix Fig. 1).
When we add the possibility of resistance gain and loss however, i.e. g>0 and l>0, the model dynamics become more complicated. Our most remarkable finding is that with this addition there may be density dependence in the outcome of competition. That is, in a closed system, one complex still eventually takes over the complete growth-space, but which one may now depend on the initial densities of the CCs (Fig. 2, compare the two panels). Specifically, a complex that starts at higher densities has an advantage, allowing it to out-compete other complexes even when these are advantaged in other ways. Such higher initial density could typically be due to earlier growth before the arrival of competitors.
The explanation for this competition effect is that higher density facilitates resistance gain (see Appendix Fig. 2). Resistant elements are transferred only if bacteria of the same CC meet, and such encounters are more likely when there are more same-CC bacteria present. If resistance is a net advantage, then the resulting higher population fraction of resistance enhances the overall fitness of the CC.
In Fig. 3, we show the outcome of competition over the range of possible starting densities of two complexes, under different conditions. As stated above, if there is no resistance transfer (g=0), the fastest growing complex always wins (A), but in case of resistance transfer we see density dependence in competitive outcome (B). If the difference in growth-rates between the CC is larger, the range of starting densities for which the slower growing CC wins is smaller (less green in Fig. 3C compared to B), and for large enough difference in growth rates, this fitness factor will always trump the population resistance benefit conferred by higher density.
The density dependent process depends strongly on the relative rates of gain and loss of resistance, for which estimates are lacking, and which might be element specific. We therefore explored different plausible combinations of these parameters in our model (see Appendix Fig. 3). The density dependence is especially strong for a medium high loss rate, or for high loss combined with high gain of the element, otherwise resistance will be either universally gained or lost for both CC populations at any prevalence level, allowing only the faster growing CC (pink in figure A3) to ever win.

Effects of continuous inflow
In a closed setting, as analysed above, eventually only one CC remains, as it will use up too much of the resource for any other CC to survive ( Figs. 2A and B, end of timeline). The hospital we aim to model is clearly not a fully closed system however, as individuals can be colonised at admittance. This can explain the long-term coexistence of complexes seen here (together with assumed CC diversity outside of the hospital); when we include inflow of different CCs in our model from elsewhere, the model predicts their co-existence (Appendix Figs. 2C and 1D, end of timeline). Yet if the inflow of MRSA from outside of the hospital is relatively little compared to the MRSA increment from growth within the hospital itself, the initial prevalence values can still determine which of the CCs dominates in the hospital subsequently, see also Fig. 3D; here, despite equal continuous introduction of both CC from elsewhere, after 10 years of competition, one of the CC Fig. 1. Example dynamics for a single CC. A: Prevalence over time of the basic MRSA resistant m (light-green) and higher resistant r-strain (dark-green) of a single CC, which has entered the hospital with low initial prevalence (m(t = 0) = 0.01 and r(t = 0) = 0.0001.) z is the density of resource available. Parameters are as at baseline for CC22 (Table 1) except inflow i = 0. B: Equilibrium prevalence of this CC as dependent on the loss-rate l of the resistant element (solid lines). The baseline loss-rate l =0.03 (used for panel A) is indicated here with a vertical dashed line. Total CC prevalence declines with increasing loss-rate, since in this setting (with high antibiotic induced death-rate, a) resistance is fitness enhancing (i.e. outweighs the cost to resistance in diminished growth, c). The equilibrium prevalence without the r-strain present is also shown (dotted lines). Note that this unstable equilibrium is lost when the mutation rate s > 0. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) .

Fig. 2.
Dynamics for two competing CCs, exemplifying density dependence. Prevalence over time is shown for the basic MRSA resistant m and higher resistant rstrains of two CCs. CC 2 has a 1% higher growth-rate than CC 1 , all other parameters are equal (and at baseline for CC22, except inflow i = 0 (see Table 1)). In A, m 1 and r 1 (light-and dark-green, together making up CC 1 ) start at double the densities of m 2 and r 2 respectively (light-and dark-purple). Panel B differs only in that these starting densities for CC 1 and CC 2 are reversed. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.). Fig. 3. Eventual outcome of competition between two MRSA CCs as determined by initial densities of both CCs. Here CC 2 (dominance in pink) has the growth advantage over CC 1 (dominance in green). For each CC, we consider only starting densities below or at the equilibrium density of this CC (as achieved without other CCs present) (hence the unequal panel sizes). Initial resistance level within each CC is assumed at equilibrium with CC density (see Methods and Appendix Fig. 2). For both CCs, parameters are as at baseline for CC22 (see Table 1) except inflow i 1 = i 2 = 0 (for panels A, B and C) and b 2 = b 1 * 1.01 (for panels A, B and D). (determined by the initial prevalences) is ~10 to 100 times as abundant as the other (note: rather than the >100 times prevalence difference after 10 years (and eventual sole survival) which occurs without external inflow). Table 2 lists all considered scenarios with their best fit to the relative prevalence data of the CCs.

Reproducing the prevalence dynamics as observed for St George's Healthcare NHS trust
With the primary model version, we can reproduce the take-over by CC22 from CC30 (Fig. 4A). In 2003, many patients at the modelled hospital became colonised with ST239 (conceivably due to the outbreak in a nearby hospital). Thus, fewer patients in St George's (and less fomites (due to any ST239 contamination from the infected), i.e. less 'resource') were left available for CC30 and CC22. Lowered prevalence of these CCs then may subsequently have impacted the level of multidrug resistance within their populations by hampering horizontal gene transfer, as explained above. When ST239 inflow dried up, CC22 more quickly recovered than CC30 did, given its higher growth rate. The relative prevalence of CC22 was further boosted by the lowered average length of stay (LOS); when the time during which others can be infected is shorter, a faster growth rate becomes relatively more advantageous. With the resulting higher relative prevalence of CC22, this complex then gained a higher population resistance level, giving it a competitive edge, which allowed it to remain dominant.
Note that in the above scenario CC22 resistance level was increased in the end as a result of its higher prevalence, which facilitates resistance transfer. Also, the good fit of this primary model version is due to the multiple stable states that are possible in this modelled system. When we remove the possibility for density dependence in competition by setting either the resistance gain rate g=0 or the loss rate l=0, this results in a quadrupling of the sum of squared differences between model and data ( Table 2). Without the density dependent process of resistance gain and loss, a reasonable fit to the data would also be possible if we make an additional assumption, namely that CC22 took an evolutionary step (Fig. 4B). If the CC22 hospital population gained a resistant element around 2004, this could also explain its growth in prevalence relative to Model output (coloured lines, for CC30 (red), CC22 (green) and ST239 (blue)) compared to the relative CC prevalences observed at St George's Healthcare NHS Trust (star points). As explained in the Methods section, the hospital system is assumed to be at steady state in 1999, meaning modelled CC levels would not change until another change occurred. In both scenarios we include two known events: an ST239 outbreak in a nearby hospital around 2004, causing a short-term high inflow of this CC, and a drop in length of hospital stay from ~6 to ~5 days in 2005. The timings of these events are indicated in the top text-bars. Panel A: Primary model version fit. Panel B: Secondary model version fit, with an additional evolutionary event assumed, causing the CC22 r-type to be introduced only in 2004, i.e. no CC22 r-type present before. See Table 2 for values of the fit parameters. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.). CC30 at that time, this additional premise in fact giving us the closest fit to the data ( Table 2). Note that in this secondary scenario, the gain of resistance by CC22 is a direct assumption rather than a model outcome.

The relative roles of ST239 and change in mean length of stay in CC22 take-over
In the primary baseline model scenario, once CC22 has reached higher prevalence and thereby higher population resistance level, it will dominate the system (Fig. 4A). The take-over by CC22 could perhaps also have been triggered by the prevalence disruption due to ST239, without the subsequent change in mean length of hospital stay. However, it probably would have taken CC22 longer to reach high levels in this case (alternative timeline shown in Fig. 5A). Conversely, the change in mean length of stay alone could have enabled relatively swift CC22 dominance take-over in the primary model (see Fig. 5B).
To explore the expected impact of these two system disruptions without the system bi-stability, we again look at scenarios without resistance transfer (g=0) (see Appendix Fig. 5). The temporary ST239 disruption could not have had a lasting impact in this case, which is most clearly shown in Appendix Fig. 5C. With a single stable balance of CC30 and CC22 prevalences, the system can only return to this state if ST239 inflow stops. The decrease in mean length of stay, although it advantages CC22 over CC30, is not big enough to cause by itself an increase in relative CC22 prevalence in the modelled hospital as large as that seen in the studied hospital (Appendix Fig. 5E). In the primary model, the system disruptions need to be amplified by the density dependent property, which emerges from the processes of CC transfer and loss of resistance. In the secondary model version, with the assumed evolutionary step in CC22 r-type inception in 2004, leaving out ST239 and/or change in average length of stay has only a minimal effect; the take-over by CC22 is here fully driven by the assumed change in CC22 alone (Appendix Fig. 5, right panels).

Fig. 5.
Counterfactual scenarios. Model output is shown for the same primary scenario fit as shown in Fig. 4A, but the model is re-run not including either one of the two known disturbance events; In panel A, the mean length of stay in hospital is kept constant (d = 0.28 also after 2005). In B, ST239 presence is not included (i 239 = 0 throughout).

Comparison with resistance data
Not only relative prevalence data, but also resistance data was documented for St George's NHS hospital. A number of isolates was available per CC per time-point, and for each isolate phenotypic resistance to eighteen antibiotics was tested ((Knight et al., 2012), summarised in Appendix Fig. 6A). Comparing these data to our model output is not straightforward however, since our model contains only two strain types per CC, rather than the many different strains with different antibiotic profiles seen in reality. We can note that in both our primary and in our secondary model versions, at baseline parameter values, we see an increase in mean CC22 resistance level after 2004, which qualitatively resembles the observed CC22 resistance data (Appendix Fig. 6B and C). Whereas in the secondary model version this increase is due directly from the assumption of late r 22 introduction, in our primary simulation it is a result of the model dynamics; higher prevalence allows for higher gain of the resistance element. In both model versions, however, CC30 loses resistance as it lowers in prevalence, which is contrary to the observation that CC30 resistance in the hospital remained stable from 2004.

What hypothetical challenger could dethrone CC22?
As a final extension of our analysis, we asked the model what type of hypothetical challenger could be expected to take over dominance from CC22 in the future. Clearly, the high density already achieved by CC22 at the hospital will make it hard for other CCs to get a foothold (see Appendix Fig. 7). Again, this is due to the density dependence in competition present in our primary model version. A successful challenger would need to be advantaged by a high growth-rate (parameter b), a considerably greater resistance to antibiotics (k), while not being limited by too high a resistance cost (c), or it would need to enter at relatively high densities.

Discussion
The aim of this work was to gain understanding of observed patterns of MRSA clonal complex (CC) dynamics, by developing and applying a new mathematical model of competing bacterial populations. Using our model, we show how theoretically density dependence may play a key role in CC competition; higher density could facilitate resistance gene build-up within the population of a CC, giving this CC a competitive edge. Furthermore, our model can reproduce the relative prevalence change over time of MRSA clones as seen actually in a single UK hospital. Here, CC22 was present from at least the end of the 1990s, but it became the dominant MRSA clone only several years later. Although datapoints are admittedly sparse, this take-over by CC22 seems to have occurred rather suddenly, going from a relative prevalence of about forty percent in 2003 to about eighty percent in 2006. CC22 was advantaged from the outset by a greater growth-rate, but it had to overcome a fitness disadvantage, compared to the earlier dominant clone CC30, in being less resistant to antibiotics. With our model we could simulate CC22 overcoming its resistance lack in two distinct ways.
First, CC22 could have become the dominant clone by having become a better competitor in a stepwise fashion (Holden et al., 2013). For example, if CC22 had gained a mutation, or if an element conferring a specific antibiotic resistance was introduced in the CC22 population around 2004 (modelled by appearance of a more resistant CC22 strain), this could have allowed CC22 to grow to dominance over the existing dominant clone CC30 at that time. However, no single new resistance type became notable in CC22 after 2003 (Knight et al., 2012). Also, take-over by CC22 cannot be explained by changed antibiotic usage in the hospital. Although an increase in mupirocin decolonisation was noted, this occurred in 2006, so after the change in dominant clone. Likewise, although nationwide increased usage of fluoroquinolones has been suspected to benefit CC22 (Holden et al., 2013;Horváth et al., 2012), resistance to fluoroquinolones was already near universal in hospital isolates from all clones in 1999 (Knight et al., 2012).
Second, we can also recreate the take-over by CC22 in this hospital without the assumption of a single evolutionary step event in CC22, or an analogous change in antibiotic pressure around 2004. In contrast to the above scenario, where resistance gain explicitly had to be assumed for CC22, increased resistance might in fact have resulted from CC22 reaching higher density. As illustrated by our model (e.g. Fig. 3), higher density facilitates resistance transfer, and this process can cause system bi-stability; the model shows how it may be that either CC30 or CC22 could become and stay the most prevalent, depending on their starting densities. In the examined UK hospital setting, the disruption by a period of high ST239 inflow from outside the hospital, from which the faster growing CC22 recovered more quickly, or the modest change in mean length of stay, which is also less problematic for a faster growing CC, could have caused the increase in CC22 abundance. Thereafter, the CC22 prevalence could increase further and achieve its resulting dominance, as its higher density then facilitated resistance gene transfer, restricted to its own CC. In other words, one of these system disruptions, or a combination of both, may have knocked the system from the observed CC30 dominant state towards the CC22 dominant state.
It seems unlikely that ST239 could have advantaged CC22 over CC30 on the longer term without this resistance transfer caused system bistability. Without bi-stability, a dominant CC would eventually return to dominance after any temporary system disruption, as this would then be the only stable state of the system. Also, the change in mean length of stay from ~6 to ~5 days seems too small to have fully caused a quick take-over by CC22, if there were no system bi-stability, which is only introduced in our model by including CC-restricted resistance loss and transfer.
As we focussed on the competition between CC, several other aspects of MRSA dynamics were treated strongly simplified in our mathematical model, or ignored. Barriers between hosts and wards were not modelled (Wang and Ruan, 2017;Bootsma et al., 2006). Within-CC-diversity in resistance was taken into account, but minimally so; the full range of actual antibiograms was simplified into a standard MRSA and a higher multidrug resistance level per CC only. In our model output, CC30 resistance levels decreased somewhat after 2004, due to density decline, contrary to seemingly stable resistance levels in the hospital. Build-up of compensatory mutations, transfer to more stable plasmids or other such concurrent evolutionary processes (Millan et al., 2014), as well as greater pressure to keep mupirocin resistance after 2006, may have kept the CC30 resistance level up. The density dependent effect could then still have occurred from the prevalence effect on CC22 alone.
Our model did not take into account how inflow of MRSA colonised patients changed over time. Including MRSA dynamics outside of hospital would have required several additional equations, with additional assumptions and parameters, for which data was lacking, such as rates of individual carriage and patient return rates (Skov and Jensen, 2009;Smith et al., 2004;Cooper et al., 2004). In fact, ultimately all connected hospitals should be modelled. Arguably, such model expansion would have been necessary to correctly model absolute prevalence rates, but our primary focus was on the process of CC competition. As a consequence, CC22 and CC30 inflow into the hospital were assumed equal and stable in our model, so that this factor did not spuriously explain the competitive dynamics between these CCs. This does not rule out the possibility that CC22 was actually overtaking other hospitals in the area, causing an increase in CC22 abundance among newly admitted patients at our case study hospital. Yet the main question would then have been shifted to how CC22 gained this advantage over CC30 in the region. Basically, we could then pose the same answers to this question as those given above for this specific hospital.
If we were to model a wider region with our model, connected hospitals (i.e. patches) with either high CC22 or high CC30 might coexist, if patient flow between them were not too high. The current model structure does not cover an explanation of why multiple CCs co-exist, as other CCs would disappear from our modelled hospital without modelled inflow from elsewhere. Heterogeneity of patches, for example differences in used antibiotics, or of patient groups among hospitals or nursing homes, likely play a role as well in keeping the full observed clonal diversity (Tilman, 1982;Miller et al., 2005;Krieger and Hill, 2018;Davies et al., 2019). Also important in this respect are the effects of random events. For example, the temporary ST239 outbreak in a nearby ICU (which could also have originated elsewhere in London) might have been due to one or a few superspreading locations, patients or devices . If it were only random events that had allowed CC22 to gain high prevalence, however, then we would still need to explain the stability of its current, long-held dominant position in the UK (Coll et al., 2017).
The simplifications in dynamics noted above gave us room to incorporate other complexities. Previous dynamic transmission models of MRSA or other antibiotic resistant bacteria have included competing strains or complexes (D'Agata et al., 2009;Kardaś-Słoma et al., 2011;Austin et al., 1999;van Kleef et al., 2017), or resistance level flux by transmissible elements (Baker et al., 2016;Freter et al., 1983). To our knowledge, this is the first MRSA model that includes both multiple competing CCs as well as resistance transfer among strains within these CCs. Importantly, due to the restriction modification system of S. aureus, transmission of genetic material rarely occurs between the different CCs (Moore and Lindsay, 2001;McCarthy et al., 2012;Goerke et al., 2004). This trait is relatively uncommon, but perhaps our model could also apply to other bacteria with clonal structure, such as Streptococcus Group B (Chen, 2019). It is this model addition which causes the density dependence in competition, since, if loss also occurs, only higher densities allow for an element to be transferred often enough for it to stay abundant within the bacterial population of a CC.
The modelled resistance transfer was kept simple by allowing for only one transmissible element per CC (i.e. two resistance levels per CC), which, together with scaling of the higher resistance level for each CC, was sufficient for our modelling purpose. Explicitly considering transfer and loss of multiple elements could render the density dependent effect even stronger. Multiple transferrable elements may have synergistic effects, and to keep each individual element abundant would require even larger population density. The population advantage might result not only from transmissible mobile genetic elements (MGEs) associated with antibiotic resistance, but also for those connected to virulence and other fitness factors (Van Wamel et al., 2006). Alternatively, or additionally, the density dependent effect we describe may play a role on different bacterial population levels. For example, on the host level, this process may play a role in explaining the fact that long-term MRSA colonisation is usually monoclonal (Cespedes et al., 2005).
The take-over of CC22 from CC30 was a nationwide phenomenon, and our theoretical model narrative for the single hospital may alternatively be considered representative for events on this larger scale. The shown importance of patient transfer networks in harmonising withinborder MRSA dynamics also suggests a country level model application (Donker et al., 2010(Donker et al., , 2017. Extrapolating to the whole of the UK, due its faster growth, CC22 was arguably more fit than CC30 already in 1999, but it took a disturbance of the system for CC22 to take over from CC30, which had the higher resistance advantage from higher prevalence at that time. ST239 seems to have been introduced to London directly from a high incidence setting in Thailand (Harris et al., 2010). Whether high inflow of this clone from Asia may have also caused disturbance in CC dynamics in other parts of the UK is unknown however, as there are few data on clonal dynamics. The decrease in mean length of hospital stay was nationwide (NHS hospital bed numbers, 2020), and alone could have triggered the switch from CC30 to CC22. Admittedly the model finding of density-dependence hinges on uncertain assumptions. However, the resulting strength of our model is that it could help explain why other CCs, such as those dominant in other countries Cockfield et al., 2007), and those noted at lower prevalence at St George's Healthcare NHS Trust (Knight et al., 2012), in the UK have not taken over since this switch to CC22 dominance (Donker et al., 2017). Nor did CC22 meanwhile take over world-wide, so it is not simply the fittest CC of all, perhaps due to its lower maximum number of resistances. According to our model, CC22 may have been stably dominant in the UK these past decades in part due to its advantageous high prevalence here.
Unfortunately, we were unable to find reference to any other longitudinal hospital dataset on clone-specific MRSA prevalence levels, which could have allowed us to better validate our model. Although we could show which assumptions are compatible with the observed phenomena, the sparseness of the data prevents us from positively affirming model correctness. Consequently, we also do not aim at any quantitative conclusions from our fits, as this would require greater certainty in parameters and model structure. How much of MRSA bacterial death otherwise induced by antibiotics is prevented by resistance genes is difficult to determine from hospital data, since it depends on the effectiveness of each type of antibiotic therapy and the frequency of use of each, but also on the average time taken to adapt treatment to specific MRSA infections. We chose to fit the resistance level per CC to the relative prevalence data to circumvent this lack of data, but our fit k parameters are a proxy for overall CC fitness differences, and should not be deemed informative on the actual effectiveness of resistance genes in avoiding antibiotic induced bacterial death.
We also unfortunately lacked data to parameterise convincingly the rates of loss and gain of transferrable MGEs in this hospital setting. As conceded, for the model results we obtained our assumptions on these parameters are crucial. For our primary scenario, we chose a loss rate which enabled the system bi-stability, and thereby density dependent competition (Appendix Fig. 3). Such a loss factor could perhaps be caused by the element regularly ending up in one daughter cell only at cell divisions (Freter et al., 1983). If the loss rate would be much smaller, all bacterial cells would soon carry the resistant element, even at low prevalence of the complex, preventing the density dependent effect. However, such consistently high resistance levels within the CCs was not seen in the hospital (Knight et al., 2012), and extensive acquisition and loss of MGEs was also observed in experimental S. aureus co-colonisation of piglets (McCarthy et al., 2014). Our demonstrated importance of resistance transfer and loss suggests then that future work should explore this heterogeneity in resistance further. For example, future work analysing resistance in multiple isolates from each colonised patient (as in Stanczak-Mrozek et al. (2015)) and categorising antibiograms of strains in hospitals over time should be used to supplement the sparse data analysis here. Supported by multidisciplinary work pairing co-culture experiments with mathematical modelling would allow for quantification of these important loss and gain rates.
Although all spread of MRSA is unwanted, understanding how policy could affect the spread of CCs differentially is relevant in that CCs might differ in their resistance potential and also their virulence (van Hal et al., 2012), and thereby in caused morbidity and mortality. Our modelling study suggest that differences between countries in main CC types present may be due to historical contingencies, and subsequent spread mostly contained within country borders by patient transfer networks (Donker et al., 2010(Donker et al., , 2017, rather than due for example only policy differences in antibiotic use. The CC that happened to be the first, locally, to incorporate an SCCmec element (Deurenberg and Stobberingh, 2008), allowing it to grow in the hospitals of that country, could claim a competitive edge due to higher prevalence from then on. This may have helped such a CC to remain dominant even when otherwise somewhat fitter competitor CCs were introduced laterour model suggests that, unless introduced at high density, only a CC with substantially higher resistance and/or growth-rate could take-over. For instance, country level antibiotic policy change might then not be expected to drive take-over by another CC more adapted to the new regime; instead, the already locally established and thereby advantaged CC could be expected to remain and subsequently adapt. However, if policy effects a stronger fitness difference between CCs, such as was observed in Hungary (Horváth et al., 2012;Conceição et al., 2007), our model does predict a switch in CC dominance, which could then affect MRSA morbidity.
In conclusion, our modelling study shows how density dependence may impact on the competition between clonal populations of MRSA, this effect potentially rendering the MRSA community in a region more stable. Thereby, instead of country level policy differences, it might be that historical contingencies mostly determine which CC has local dominance.
No funders played any role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Declarations of Interest
None.