Mathematical Modelling of Pseudomonas aeruginosa L-forms Reveals Complex Interplay Between Host Defence Mechanisms and Putative Treatments

The bacterium Pseudomonas aeruginosa has been shown to undergo a morphological transition akin to L-forms under exposure to antibiotics, a process which may contribute to persistent infections. With the further consideration of antibiotic-resistance mechanisms, this transition renders the design of effective treatment strategies challenging. Through a mathematical model, we illustrate that additionally incorporating the complexities of the host immune response can render somewhat surprising predictions from the simulations. In particular, scenarios arise whereby the addition of a treatment strategy to directly target the L-forms results in a worsened infection, while in others this treatment could turn an antibiotic-resistant infection from persistent to treatable. The study highlights the importance of understanding the in vivo interplay between immune cells and pathogens for successful treatment design.


INTRODUCTION
Pseudomonas aeruginosa is a Gram-negative pathogenic bacteria, categorised as critical by the World Health Organisation in the search for novel therapeutics to combat antimicrobial resistance (WHO, 2017). The β-lactam antibiotic meropenem works by targeting the cell wall of P. aeruginosa and inhibits its cell wall synthesis by binding to penicillin binding receptors. This mechanism of action results in a breakdown of the cell wall, which leads to filamentation and the formation of spheroplasts, which eventually lyse destroying the bacteria (Heinemann et al., 1998). However, in the presence of meropenem, P. aeruginosa has been observed to undergo a rapid morphological change: it changes from its usual rod-shaped cell to a spherical-shaped cell with a defective cell wall (Monahan et al., 2014). This sub-population of cells has the capacity to evade the effects of meropenem through the absence of crucial targets in this cell wall. These spherical-shaped cells remain dormant in vitro until the removal of the antibiotic triggers their transition back to proliferating rod-shaped cells, thus drawing analogies with persister cells and L-forms that have been implicated in chronic infections (Errington et al., 2016).
Persistence differs from genotypic drug resistance as, unlike the latter, it is not inheritable and is thought to be the result of an intrinsic survival mechanism. Additionally, unlike genetically resistant cells, when a persister cell reverts back to its original cell-type it is once again susceptible to antibiotics (thus, this is often referred to as antibiotic-tolerance). The spherical-shaped cells into which P. aeruginosa can transform are not unlike persister cells but the potential fragility of the P. aeruginosa spherical cells may lead to a high natural death rate, and perhaps as a result of this they are not always referred to as genuine persister cells. It is nonetheless quite plausible that they will play a role in chronic infections (Monahan et al., 2014). Thus while in this study we do refer to these spherical cells as L-forms/ persister cells for ease of comparison to similar bacterial phenotypes, this potential distinction should be noted.
It has been postulated that antimicrobial peptides targeting the spherical cells could be a useful adjuvant to antibiotics to curtail chronic infections, and we previously developed a mathematical model that illustrated the feasibility of this in an in vitro setting (Spalding et al., 2018), agreeing with the experimental results of (Monahan et al., 2014). In the current study, we extend this model to consider the putative effects of the immune response on such a treatment, including cases where the bacteria may be additionally genetically resistant to antibiotics. The crucial interplay between the host defense and pathogen population results in often counter-intuitive outcomes. While there are scenarios where the immune system may be capable of clearing the persistent spherical cells without additional treatment (similar to the findings of (Gjini and Brito, 2016) where the immune response was found to potentially be critical in clearing antimicrobial resistant cells), there are also biologically plausible scenarios where adding AMPs could surprisingly have an adverse effect on the ability of the immune system to respond and therefore on the treatment outcome. For optimal treatment designs, this study highlights the need for mathematical models to consider the interactions of the host defense mechanisms with the pathogen and treatment load with a great deal of care.

MODEL FORMULATION
The model describes the growth dynamics of a multi-strain P. aeruginosa population with the influences of antibiotics and the immune system at time t. The following subpopulations are used with i = S, R denoting antibiotic-susceptible and -resistant variants respectively: • Antibiotic-free rod-shaped bacteria, H i (t); • Rod-shaped bacteria with internalised antibiotic, I i (t); • Spherical (L-form) cells, L i (t); with extracellular antibiotic (A) included as a parameter, and the immune response captured through a variable, P(t), for phagocytes.
We assume the bacteria grow according to logistic growth, (with specific growth rate r and carrying capacity k), and assign a fitness cost parameter, c, to the growth of a resistant strain, as is done, for example, in (Bergstrom et al., 2000;Ternent et al., 2015;Baker et al., 2016). We neglect natural death for all subpopulations bar spherical cells (L i ) that are likely more prone to this fate due to their fragility (occurring at rate ψ). Figure 1 illustrates how the subpopulations respond to the presence of antibiotic. Cells are rod-shaped in their 'natural' state.
In the presence of antibiotic, rod-shaped cells either.
• internalise the antibiotic (which prevents cell division and can potentially kill the cells), or • transform to spherical L-forms (cell division is also obstructed here, see for example (Mercier et al., 2014;Cross et al., 2019;Chikada et al., 2021) where proliferation either does not occur or is seen to only occur in very specific environmental conditions).
Antibiotic-internalisation is modelled using a non-linear saturating term with maximum transition rate, ], and antibiotic concentration needed for half-maximal internalisation, A 50 . To represent resistance, we model efflux pump upregulation as simply a decrease in the internalisation rate, ], and define ω ∈ (0, 1] to be the efflux pump upregulation parameter (resistance may be partial; w = 1 signifies maximal resistance).
Rod-shaped bacteria with internalised antibiotic may lyse due to antibiotic effects as the antibiotic inhibits cell wall synthesis. We define ρ to be the death rate of rod-shaped cells due to antibiotic effects and assume that this death rate is the same for both strains since resistance is instead captured via a decrease in internalisation. We note that A (antibiotic concentration) is not needed explicitly in these death terms as it is implicitly included in the variables I S and I R .  Table 1. Here, D represents dead cells. Spherical cells can lyse naturally due to fragility and bacteria that have internalised antibiotic can die due to the effects of the meropenem.
Following the microscopy images in Figure 2, we assume that both strains are able to transition to spherical cells following exposure to the antibiotic. The transition is modelled using a saturating nonlinear term with maximum transition rate, γ, and concentration needed for half-maximal transition, T 50 .
For simplicity, we only consider continuous administration of the antibiotic in this study and hence the cells cannot transition back to the rod-shape (see (Spalding, 2018) for alternative administration strategies).
To capture an immune response, we assume that all bacterial subpopulations are susceptible to phagocytosis. The immunogenic properties of L-forms remain relatively unknown but it has been postulated that they will differ to those of actively growing cells (Domingue and Woody, 1997;Onwuamaegbu et al., 2005;Errington et al., 2016) where the cell envelope is highly immunogenic: both the outer membrane and peptidoglycan layer have constituents that act as pathogenassociated molecular patterns (PAMPs). These molecules are recognisable to the immune system and phagocytes will respond accordingly. Therefore, since the cell wall is disrupted in the spherical cells it is not clear whether or not the spherical cells will incur the same immune response as cells that harbour an intact cell wall. In line with other mathematical models of phagocyte recruitment (Smith et al., 2011;Ternent et al., 2015), we model phagocyte population growth using logistic growth, and assume that the morphology of a bacterium could affect its immunogenic properties. We define the maximum recruitment rate κ i and phagocytosis rate ϕ i , where i H, L { }, with subscripts H and L referring to the rod-shaped and spherical bacteria, respectively. The carrying capacity of phagocytes is defined as P MAX .
P. aeruginosa bacteria can produce toxins that induce apoptosis of neutrophils and other immune cells (Usher et al., 2002). We model bacteria-induced apoptosis using a linear density-dependent term with rate, μ. For simplicity, we assume that bacteria-induced apoptosis occurs at the same rate for all bacterial subpopulations. Additionally, we define a clearance rate of phagocytes, ξ.
We consider the introduction of resistance into the environment via spontaneous mutation or cross contamination, i.e. a resistant population is introduced at some point (T R ) over the course of the simulated treatment. If we are considering a resistant population that has arisen due to mutation then our initial population size (H R0 ) would be small, whereas a population that appears from cross contamination could be much larger. Our aim is to understand the conditions under which a resistant population would grow and to investigate how treatments could work in conjunction with the host response to control the emergence of resistance and clear the bacterial infection.
FIGURE 2 | Cells of the P. aeruginosa (A,B) meropenem-susceptible PA1008 strain and (C,D) partially resistant PA1004 strain rapidly convert to a spherical shape after adding (A,C) 0.5 μg ml −1 and (B,D) 2 μg ml −1 of meropenem at 0 h before reverting back to a bacillary form. Left panels: transmitted light, right panels: fluorescence microscopy. In fluorescence images, green and red colouring indicates viable and lysed cells respectively (see (Spalding et al., 2018) for details).
Frontiers in Systems Biology | www.frontiersin.org June 2022 | Volume 2 | Article 899990 Since we are exploring the limiting case of constant antibiotic administration, antibiotic is included as a parameter in the resulting model: with initial conditions, The definitions of all variables and parameters are displayed in Tables 1 and 2 respectively. A default parameter set was chosen using suitable literature, the parametrisation results from (Spalding et al., 2018), and our numerical investigations to ensure the model behaves realistically in the absence of an immune response or drugs (i.e., the infection takes hold in the absence of either, and the immune response alone cannot clear the infection). The default parameter values are shown in Table 2, although a range of values is used throughout and additional parameter sets have been tested to account for the inevitable uncertainty in the default choice.
To produce the simulations in this paper, as in (Ternent et al., 2015), we include the assumption that bacterial growth will cease when the population reaches a low critical threshold, i.e.  environmental pressure and the stochasticity of biological processes would cause the rod-shaped population to die out if the population decreased to a sufficiently low level. When solving the model we include the assumption that growth will stop once the rod-shaped population reaches 50 cells cm −3 and further assume that the rodshaped bacteria will no longer transition to the spherical form.
In the following results, we begin by considering the case where spherical cells and rod-shaped cells have identical immunogenic properties (κ H = κ L , ϕ H = ϕ L ) before exploring the more complex consequences of these rates diverging by cell type.

Assumption: Rods and Spheres Have Identical Immunogenic Properties
We first consider the scenario whereby immune cells do not distinguish between rod-shaped and spherical cells.

No Treatment
To illustrate our parameter choice for the immune-related parameters, initially we consider bacterial growth of a fully susceptible strain, prior to treatment (A = 0), such that the initial conditions are Eq. 8 with, H S0 10 4 , H R0 0.
(9) Figure 3A displays the solutions of the model using the default parameter set with A = 0. All bacteria remain in the rodshaped form since there is no antibiotic present to induce changes. κ H (recruitment rate) and ϕ H (phagocytosis rate) are chosen such that the immune response cannot clear the infection alone.
We now consider a single strain that is partially-resistant to the antibiotic and solve the model using initial conditions Eq. 8 and In (E), H R0 = 100 and the simulation is run for 5,000 days (as the steady state is reached on a unrealistically long timescale). In (F), the resistant bacteria have no fitness cost (c = 0) to consider the limit case, and steady state results (R ∞ are shown).
Frontiers in Systems Biology | www.frontiersin.org June 2022 | Volume 2 | Article 899990 to reflect the population being all antibiotic-resistant rod-shaped bacteria ( Figure 3B). Although the initial population size is the same, the total bacterial load is lower than the susceptible case due to the fitness cost, c, of resistance. In Figures 3C-F, we illustrate how the initial number (H R0 ) and time at which resistance first appears (T R ) affects the long-term resistant bacterial load in a mixed susceptible-resistant infection (H S (0) = 10 4 and H R0 is varied). The resistant population inevitably dies out eventually except in certain cases where the fitness cost is very low (e.g., lower right section of Figures 3C,D). Fitness cost values of the size needed for extinction in Figure 3 are extremely small in comparison to the interval of realistic values provided in (Subbiah et al., 2011;Vogwill et al., 2016) and although these values are related to plasmid-mediated resistance, we can use them as an indication of the potential size of a fitness cost, i.e., it is likely such small values of c are unrealstic. Nevertheless, we include them in the sensitivity analysis for completeness. Figures 3C,D indicate that the initial condition will determine the size of the resistant population. In contrast, Figures 3E,F indicates that T R does not significantly affect the population size of the resistant population in the long term when compared to the effects of the fitness cost and initial condition. Results for the total number of susceptible bacteria are not shown as they are not significantly affected, regardless of the characteristics of the resistant population.

Antibiotic Treatment
We first consider a population of only susceptible bacteria (initial conditions Eq. 9)-see green dashed lines in Figures 4A-D.
Within the first few days we see a rapid decrease in rod-shaped bacteria with no internalised antibiotic (H S ) as the antibiotic binds within the bacteria or they shed the cell wall. The presence of the bacterial population initiates an immune response (P). As the antibiotic treatment continues we see all subpopulations decrease until the infection is cleared, albeit after an undesirably lengthy treatment period. We note that our results in (Spalding et al., 2018) suggested that the morphological transition could enable an antibiotic-susceptible population to survive a constant dose of antibiotic; the results in Figures 4A-D suggest that the additional effects of the immune response could aid the antibiotic and clear an antibiotic-susceptible infection. This result can be directly related to the study by (Ankomah and Levin, 2014) in which they suggest that the presence of a subpopulation with phenotypic resistance (e.g., persister cells) would not cause treatment failure in cases where the rate of phagocytosis is sufficiently high. These results highlight the importance of acknowledging hosts defences when considering mechanisms of antibiotic evasion. Conversely, if the population is entirely antibiotic resistant (initial condition Eq. 10 and blue dot-dashed line in Figures  4A-D), the resistant bacteria would reach a non-zero steady state, even after continued exposure to the antibiotic and in the presence of an intact immune response.
Using the mixed initial conditions, H S0 10 4 , H R0 100, Figure 4E displays the steady state values for the resistant population when the time of resistance introduction is varied against the initial number of resistant bacteria-unlike the antibiotic-free case in Section 3.1.1, neither variable affects the steady state population. In fact, neither T R nor R 0 has any impact on the steady state population level when compared with the fitness cost, c, or the efflux parameter, ω (see Figures 4F-I). Figure 4J displays the results if we vary the fitness cost, c, and efflux upregulation parameter, ω. Two clear regions of behaviour emerge for low and high values of c. When we have a low fitness cost (c < 10 -2 for our parameter set), we find that changes to the efflux parameter have a large impact on the steady state of resistant bacteria whilst the fitness cost has little effect. This is reversed for higher fitness costs that prevent the long-term survival of the resistant population.
Note that Figures 4G,I both suggest that the threshold fitness cost, below which resistant bacteria would survive, is higher in this case than in the antibiotic-free environment ( Figures 3C-E). This arises due to the niche created by the mass killing of the susceptible bacteria by the antibiotic. These results agree with the concept that antibiotic exposure causes a selective pressure that could enable resistant populations to emerge that would not survive in an antibiotic-free system.

Antimicrobial Peptide Treatment as an Adjuvant to Antibiotics
The steady-state survival of bacteria in the above section is predicted to be due to efflux-mediated resistance rather than the morphological change, since the immune response ultimately tackles the L-forms. This ability to create L-forms, however, does enable susceptible bacteria to survive for lengthy periods of time before being cleared. Lengthy treatment periods would be costly and provide an ideal environment for further resistance mechanisms to develop and leave the patient vulnerable to more threatening infections. In (Monahan et al., 2014) antimicrobial peptides were suggested as a suitable therapy to target L-forms and our results in (Spalding et al., 2018) supported the use of antibiotics and AMPs when the bacteria was in isolation with the antibiotic (i.e., in the absence of an immune system).
We simulate the addition of AMPs through an increase in the lysis rate of spherical cells in Figure 5 for an entirely antibioticsusceptible infection. If the antibiotic is administered alone (red dashed lines) then we predict that the population comprises only spherical cells towards the end of the treatment. Therefore, the antibiotic at the later stages is not bactericidal but instead preventing the spherical cells from reverting back to rods. The remaining spherical cells are then eradicated through natural clearance and phagocytosis. In theory, the AMPs should have the ability to clear this population of spherical cells faster. However, their potency has a detrimental effect on the immune response: if spherical cells are killed by the AMPs too quickly, the early immune response reduces accordingly and rod-shaped cells survive for longer. Eventually, the immune response recovers and the rod-shaped population is cleared, but the treatment duration is still lengthy (note that the sharp dropoff in population densities is due to the threshold below which cells can no longer proliferate-see end of Section 2). Thus in essence, the addition of AMPs really only affects the composition of the population: for higher values of ψ (addition of more peptides), the population towards the end of the treatment comprises mostly rod-shaped cells with no internalised antibiotic; for lower values the spherical cells are the surviving population. We recall that due to the compromised cell membrane the spherical cells could be less virulent than the rod-shaped cells. If so, using AMPs could result in a more threatening infection. Figure 6 depicts the minimum time needed to eradicate the population varying ψ. We estimate the optimal death rate of the spherical cells to be ψ ≈ 0.25 day −1 for our default parameter set. Beyond this synergistic potency level, the AMPs could become ineffective and possibly detrimental to treatment efforts.
If we simulate a combined dose of antibiotics and AMPs as a treatment for a solely antibiotic-resistant population then the results, detailed in Figure 7, suggest similar findings to those obtained for the susceptible strain. Antibiotics alone cannot clear the infection and as we increase the spherical cells' death rate from the default parameter value we predict fewer spherical cells and a larger rod-shaped bacteria population. Overall, the steady state of the total population is similar, regardless of the efficacy of the AMPs and their influence on the immune response. As before, it is the composition of the population that is mostly affected, not the total size.

Assumption: Rods and Spheres Are Phagocytosed at Different Rates
As mentioned in Section 2, it is unknown how well immune cells phagocytose spherical cells. Phagocytosis may occur at a faster rate since the cells have a depleted cell wall. However, the depletions may incur loss of targets for the immune cells and therefore phagocytosis could occur at a slower rate. Investigation into higher phagocytosis rates for the spherical cells (ϕ L > ϕ H ) indicated that the bacteria would be cleared faster either with antibiotics alone or with AMPs as an adjuvant to antibiotics, than in the case where ϕ H = ϕ L (results omitted for brevity). This makes sense biologically since the 'persistence' mechanism is being targeted at a higher rate by the immune system and there is nothing interfering with the immune response. We turn, instead, to the case where ϕ L < ϕ H , i.e. the immune cells phagocytose spheres at a slower rate than they do rods.

Antibiotic-Susceptible Infection
Decreasing the spherical cell phagocytosis rate, ϕ L , leads to the results in Figure 8. The model predicts larger populations of spherical cells as the phagocytosis rate is decreased from its default value. The larger spherical population causes the total population of bacteria to increase and this initiates a stronger immune response as more phagocytes are recruited. Additionally, since fewer spherical cells are phagocytosed, fewer phagocytes undergo apoptosis. Consequently, and perhaps counter intuitively, the higher number of phagocytes causes the rodshaped populations to be cleared faster when the phagocytosis rate of the spherical cells decreases. The spherical cells are susceptible to natural death and replenishment depends on the size of the rod-shaped population. Hence we see a decline in the spherical population once the rod-shaped population has depleted.  Table 2 and varying the spherical cell death rate, ψ, from the default value of 0.017 day −1 . An increase in ψ represents the addition of antimicrobial peptides.  Table 2 and varying the spherical cell death rate, ψ, from the default value of 0.017 day −1 . An increase in ψ represents the addition of antimicrobial peptides.
Frontiers in Systems Biology | www.frontiersin.org June 2022 | Volume 2 | Article 899990 9 The qualitative result for the total population remains the same regardless of the phagocytosis rate: we predict total bacterial clearance, even if we assume that the spherical cells are immune to phagocytosis (ϕ L = 0 cyan solid lines). The treatment duration needed to clear the infection increases as the phagocytosis rate decreases, and we could expect a higher bacterial load over the course of the treatment, albeit with the population dominated by spherical cells (note that antibiotics are required over the full course to prevent the spherical cells transitioning back to rod-shaped cells). Similarly, the addition of AMPs to target the spherical cells from the start of treatment encounters the same issues as seen in Section 3.1.3: the ensuing early reduction in the immune response can actually lengthen the treatment duration (details can be found in (Spalding, 2018)).

Antibiotic-Resistant Infection
Figure 9 displays the model solutions for an antibiotic-resistant infection if we vary the phagocytosis rate of the spherical cells, ϕ L , assuming that ϕ L < ϕ H . The results suggest similar findings to the susceptible strain: as the spherical cell phagocytosis rate decreases from the default value, the spherical population increases and initiates a stronger immune response. However, in the case of an antibiotic-resistant infection, the increase in phagocytes could kill off the rod-shaped population that would otherwise survive if the immunogenic properties were identical for all cell types (compare with the blue dot-dashed lines of Figures 4A-D). Therefore, we predict that a resistant strain infection could become treatable with antibiotics if the spherical cell phagocytosis rate is slower than that of the rod-shaped cells.
The results in Figure 9 are counter-intuitive and it is only through consideration of the immune response that we are able to obtain these predictions. The results suggest that it may not be beneficial for the resistant bacteria to activate the morphological transition at a site of infection if the spherical cells are phagocytosed more slowly than the rod-shaped cells. If the resistant population suffers as a result of transitioning to the spherical form then it is possible that they may evolve to not utilise this mechanism; we will consider this in the discussion.
We next assume that the spherical cell phagocytosis rate is lower than that for rods and simulate the use of AMPs along with the antibiotic. If sufficiently potent the AMPs initially cause a decrease in spherical cells and therefore there are fewer phagocytes recruited than when the antibiotic was added alone ( Figure 10). Consequently, the rod-shaped population is not overwhelmed by the immune response and the population survives treatment. These results suggest that AMPs could provide a more favourable environment for the rod-shaped population and again prove detrimental to treatment success. Figure 11 illustrates the results when we vary the phagocytosis rate for the spherical cells in a mixed strain infection. For both strains, the results indicate similar findings to those obtained when considering single strain populations. As the spherical cell phagocytosis rate is decreased, a larger spherical population that comprises mostly antibiotic-susceptible spherical cells emerges. The increased spherical population initiates stronger immune cell recruitment and subsequently the rod-shaped subpopulations die out faster. Eventually the susceptible population is cleared after continued antibiotic usage. As previously illusrtated in FIGURE 9 | The variable solutions to system Eqs. 1-7, with initial conditions Eq. 10 (entirely resistant infection), default parameter values shown in Table 2 and varying the spherical cell phagocytosis rate, ϕ L , from the default value of 2.4 × 10 -4 (cells cm −3 ) −1 day −1 .

Mixed Strain
Frontiers in Systems Biology | www.frontiersin.org June 2022 | Volume 2 | Article 899990 Figure 9, the rapid morphological transition could also potentially suppress the emergence of resistance through this boost to the phagocyte population. However, we still estimate a lengthy treatment period and conclude that a successful treatment would need to clear the spherical populations to prevent either strain recovering after the treatment is stopped. If we simulate the use of AMPs along with the antibiotic, we obtain the results in Figure 12. The effects of the AMPs cause a decrease in susceptible spherical cells and therefore there are fewer phagocytes recruited. When only antibiotics are used, the large phagocyte population (recruited due to the large susceptible spherical population) showed potential to delay the emergence of the resistant population. By adding AMPs we are providing a more favourable environment for the resistant strain (phagocytes are the main means to kill these cells) and we see the resistant rod-shaped cell population flourish for sufficiently large ψ. The increase in spherical cell death should negatively impact both of the spherical populations but the rod-shaped resistant population becomes so large that the AMPs do not prevent growth of the resistant population. Our predictions indicate that in the case of a multistrain infection with a spherical population that is less susceptible to phagocytosis, the use of AMPs could be detrimental and may not pose a suitable treatment in combination with antibiotics.

Assumption: Rods and Spheres Recruit Immune Cells at Different Rates
As the outer membrane may still be intact in the spherical cells, it may be reasonable to assume that shedding the cell wall is more likely to affect the recruitment rate than the phagocytosis rate. PAMPs are often found in the peptidoglycan layer of the cell membrane and the outer membrane and therefore it is possible that the spherical cells will have fewer PAMPs and be less immune stimulatory (κ L < κ H ).

Antibiotic-Susceptible Infection
With κ L < κ H , in the initial stages of treatment, the lower recruitment rate of the spherical cells causes fewer phagocytes to be recruited to the infection ( Figure 13). As a result of the smaller phagocyte population, all bacterial subpopulations increase as the recruitment rate of the spherical cells decreases. If the spherical cell recruitment rate decreases sufficiently, treatment now fails.
For Figure 14 we adopt the phagocyte recruitment rate of the spherical cells to be κ L = 1 × 10 -3 (cells cm −3 ) −1 day −1 (i.e., treatment fails) and investigate whether the infection can be cleared using a combination therapy of antibiotic with AMPs. The results suggest that the addition of AMPs could lead to fewer spherical cells. However, as the AMPs only target the spherical cells and the spherical cells have a low phagocyte recruitment rate, we find that the remaining subpopulations are unaffected by the addition of AMPs (indeed the lines are indistinguishable). Thus while, adding AMPs could potentially help clear the infection, this could still require a lengthy treatment duration.

Antibiotic-Resistant Infection
The consequences of the spherical cells having a lower phagocyte recruitment rate than the rod-shaped cells cause identical changes FIGURE 10 | The variable solutions to system Eqs. 1-7, with initial conditions Eq. 10, default parameter values shown in Table 2, with ϕ L = 1 × 10 -5 (cells cm −3 ) −1 day −1 and varying the spherical cell death rate, ψ, from the default value of 1.7 × 10 -2 day −1 . An increase in ψ represents the addition of antimicrobial peptides.
on an antibiotic-resistant infection to the results described in Section 3.3.1 for the antibiotic-susceptible case. For brevity we do not include the corresponding simulations.

Mixed Strain
In a mixed strain infection (Figure 15), the model predicts similar dynamics to the single strain up to around day 20: all susceptible subpopulations grow to higher levels as we decrease spherical-cellinduced phagocyte recruitment. Due to the lower recruitment rate, we initially predict a weaker immune response and as a result we see faster growth of the resistant population. However, as the resistant population grows, the immune response increases. Consequently, the antibiotic and increased immune response can clear the susceptible population that previously persisted in the single strain simulations. This predicted behaviour suggests that if recruitment is lower for the spherical cells then we could expect the presence of a resistant strain to impact the survival of the susceptible strain negatively whilst promoting its own growth.
If we simulate the combination of AMPs and antibiotics then we obtain the predictions in Figure 16. As expected, the AMPs are successful in clearing the spherical cells, however, they do not affect the rod-shaped cell populations and therefore we could still expect a dominant resistant population to emerge.

DISCUSSION
We have formulated a model that describes the growth dynamics of P. aeruginosa at the site of infection with the inclusion of a morphological transition, immune response, antibiotics and AMPs. We summarise our results for all treatment combinations and all explored immunogenic characteristic settings in Table 3.
If we assume that the immunogenic properties are identical between cell types then the results suggest that antibiotic treatment could clear a P. aeruginosa infection where L-forms can emerge only if it was purely susceptible to the antibiotic. In (Spalding et al.,FIGURE 11 | Numerical solutions of Eqs. 1-7, with initial conditions Eq. 11 (mixed infection), default parameter values shown in Table 2 and varying the spherical cell phagocytosis rate, ϕ L , from the default value of 2.4 × 10 -4 (cells cm −3 ) −1 day −1 .
Frontiers in Systems Biology | www.frontiersin.org June 2022 | Volume 2 | Article 899990 2018) we predicted that a continuous dose of antibiotic would not clear a persistent infection and after including an immune response within the model we see that treatment success is now predicted. This result supports the conclusions in (Ankomah and Levin, 2014) and (Gjini and Brito, 2016) and emphasises the importance of including the immune response when considering persistent populations, as host defences will play an important role in the treatment outcome. If a resistant strain is present then we predict that the antibiotic treatment would fail, as expected. The varied results in Table 3 indicate that the immunogenic properties of the spherical population could affect the survival of an infection and change the treatment outcome. If we examine the predictions for bacterial growth during antibiotic treatment and vary the immunogenic properties then we recognise two interesting scenarios.
Firstly, if the recruitment rate was slower for the spherical cells, we predicted that the susceptible population could be untreatable with antibiotics alone, and additional treatment types would be needed for infection eradication.
Secondly, we found that if the phagocytosis rate of the spherical cells is lower than the rod-shaped cells, a resistant infection could be treatable with antibiotics alone. The resistant bacteria cause a large influx of phagocytes that enable total bacterial clearance. This compelling result seems counter-intuitive and leads us to question whether the resistant bacteria would utilise the morphological transition in vivo if it was detrimental to the overall population. It is possible that the resistant bacteria could adapt and remain as bacillary cells if the selective pressures of the immune response caused by the antibiotic-induced morphological transition negatively impacted bacterial fitness.
Our results suggest that a susceptible strain could be cleared using antibiotic alone but that a lengthy treatment period would be necessary. After considering dual treatments we can formulate predictions about the clinical potential for AMPs. In (Spalding et al., 2018) we produced results that strongly supported the combined use of meropenem and AMPs to clear an antibioticsusceptible bacterial population in vitro. Whilst we do not have FIGURE 12 | The variable solutions to system Eqs. 1-7, with initial conditions Eq. 11 (mixed infection), default parameter values shown in Table 2, with ϕ L = 1 × 10 -5 < ϕ H (cells cm −3 ) −1 day −1 and varying the spherical cell death rate, ψ, from the default value of 1.7 × 10 -2 day −1 to mimic the addition of AMPs. experimental data to indicate the analogous result for a resistant strain, it is feasible that a resistant population could also be cleared with the antibiotic and AMPs in isolation. However, after including host defences, we find that the results do not strongly support the use of AMPs when used as part of a combined continuous therapeutic strategy. Our results suggest that the success of a combined treatment of an antibiotic and AMPs would depend on the effectiveness of the AMPs and the composition of the infection and, in many cases, the addition of AMPs was predicted to have a detrimental effect on the immune response and hence worsen the infection. This opens up an important question of timing-if AMPs were to be added only after the rod-shaped bateria had been cleared, this sequential treatment would ensure treatment success. Mathematical models will surely play a role in determining such personalised strategies.
There are, however, simulated scenarios where AMPs were predicted to improve treatment, and this work is in no way intended to condemn their use. Rather, similar to (Gjini and Brito, 2016), we aim to highlight the need for a better understanding of the interplay between the host, the pathogen and the treatment to design therapeutic strategies (Wheatley et al., 2021). Additional treatment types, such as phage therapy or anti-virulence drugs could also be considered for combination therapies (Spalding, 2018;Theuretzbacher and Piddock, 2019), and both experimental and theoretical studies have emphasised the need to consider careful timing of treatments and potential antagonism ( (Ternent et al., 2015;Cogan et al., 2016;Gjini and Brito, 2016;Rezzoagli et al., 2020) for example).

SUMMARY
We have presented a model of P. aeruginosa and its responses to meropenem and host defences, emphasising the sometimes FIGURE 15 | The variable solutions to system Eqs. 1-7, with initial conditions Eq. 11, default parameter values shown in Table 2 and varying the spherical cell phagocyte recruitment rate, κ L , from the default value of 0.5 day −1 .
surprising outcomes. While the model would of course benefit from experimental data against which we could validate and narrow down realistic scenarios, we have clearly demonstrated the importance of considering interactions with the host immune response when designing effective treatment strategies to counter chronic infections.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
CS, SJ and AM-K conceived and designed the analysis; EK collected the data; CS and SS contributed data or analysis tools; CS performed the analysis; CS, SS and SJ wrote the paper.
FIGURE 16 | The variable solutions to system Eqs. 1-7, with initial conditions Eq. 11, default parameter values shown in Table 2, with κ L = 1 × 10 -3 < κ H day −1 and increasing the spherical cell death rate, ψ, from the default value of 1.7 × 10 -2 day −1 to mimic the addition of antimicrobial peptides.