A within-host infection model to explore tolerance and resistance

One purpose of our work was to develop a mathematical model capable of reproducing the wide variety of WHD documented to date. Our model aims at being general, but we used the recent empirical work of Duneau et al., 2017a on Drosophila melanogaster to facilitate explanations and interpretations. This previous work documents chronic and acute diseases, some benign, while others are lethal. But most importantly, it also depicts situations where an infection initiated in carefully controlled conditions can bifurcate, with some of the hosts developing a form of chronic infection when others die within a few days (as illustrated in Figure 1).

Such bifurcations in dynamical systems are most often by-products of bistability. The few available examples of bistable WHD models in the literature are all modifications of the Lotka-Volterra predator-prey model (Pujol et al., 2009; Souto-Maior et al., 2018). Mayer et al., 1995 demonstrated that in such models non-linear functional responses for both pathogen dynamics and host defense regulation could yield bifurcations between chronic infections with high and low equilibrium pathogen load.

Other models have proven that the non-linearities that create bifurcations could originate from feedback between pathogen proliferation and immune defense production. van Leeuwen et al., 2019, for example, demonstrated that a model that explicitly describes the amount of resource intestinal pathogens diverted from their host can be bistable, a result later confirmed by Yu et al., 2021. Ellner et al., 2021 showed that situations where pathogens actively destroy immune defenses and also create bistability in disease outcome. All these models have in common that they predict bistability where some infections are cleared (or almost cleared, in the case of Ellner et al., 2021 conceptual model). In that sense, they do not reproduce the situations documented by Duneau et al., 2017a, where some infections remain chronic while others are lethal.

In our research, we amalgamated Mayer et al., 1995’s concept of nonlinear modulation in defense production with the hypothesis that pathogens cause damage that can hinder the host’s ability to combat infection. We customized our model to not only forecast pathogen dynamics, but also estimate the host’s survival prospects.

Let ${\displaystyle x}$ be the pathogen load (e.g. bacteria), ${\displaystyle y}$ the level of defense (e.g. anti-microbial peptides) produced by the host immune system, and ${\displaystyle z}$ the damage caused by the infection (see Figure 2A). Whether the host suffers from severe infection, chronic yet mild infection, or clears pathogens, depends on its genetically determined capacity to fight or tolerate pathogens, but also on the inoculum size and on the host physiological state at the start of the infection. In our model, the inoculum size is the initial value of the pathogen load (${\displaystyle x_0}$). The physiological state of the host at the start of the infection is described by the initial level of defense (${\displaystyle y_0}$) and of damage (${\displaystyle z_0}$). The initial level of defense can be set by constitutive immunity, when a healthy host is at its homeostatic state (i.e. ${\displaystyle y_0={\tilde{y}}_h}$). It can also be lowered, when the host is weakened by environmental challenges. Similarly, before the infection has started, the host suffers from an initial level of damage which origin is not the studied infection. This level is ${\displaystyle {\tilde{z}}_h}$ when the host is at its homeostatic state, but can be increased in case of environmental challenges (e.g. infection starting from a wound).

Figure 2

Download asset Open asset

A set of parameters (described in Figure 2A) characterizes the genetic interaction between a host and a pathogen. These parameters reflect both the genetically determined capacity of the host to fight an infection and the genetically determined aggressiveness of the pathogen it is fighting against. These parameters govern the dynamics of ${\displaystyle (x,y,z)}$ over time ${\displaystyle t}$, as described by the following differential equations:

For the sake of simplicity, the equations of system (Equation 1) are written dimensionless, ${\displaystyle x}$ being a fraction of the pathogen’s carrying capacity within the host and time units corresponding to the pathogen’s generation time. The equation governing the dynamics of pathogens is that of the evolution of a prey population in a classic Lotka-Volterra model (Lotka, 1923; Volterra, 1928). The pathogens are destroyed by the defense at a rate ${\displaystyle \delta y}$, with ${\displaystyle \delta }$ representing the efficacy of immune molecules or cells at killing pathogens.

We considered that damage should increase with pathogen load, but that it might also increase with the level of defense that is produced (i.e. immunopathology, González-González and Wayne, 2020; Graham et al., 2005; Rommelaere et al., 2024). We assumed that damage are caused by pathogens at a rate ${\displaystyle \omega }$ and by defense at a rate ${\displaystyle \eta }$. We further reasoned that damage are repaired at a constant rate ${\displaystyle \xi }$. It is now well established that repairing mechanisms are in fact modulated during the course of an infection (Martins et al., 2019; Pucillo and Vitale, 2020). For simplicity, and because the details of how this regulation responds to damage and defense production are still unknown, we deliberately neglected this modulation. We also made the simplifying assumption that damage caused by pathogens and by immune responses is of the same nature, allowing it to be repaired at the same rate. Finally, we postulated that the host cannot sustain infinite damage, resulting in death if ${\displaystyle z}$ exceeds the threshold value ${\displaystyle z_d}$.

We assumed that the time variation in the level of defense is set by the balance between the production of new molecules or cells (${\displaystyle F(x)G(y,z)}$) and the spontaneous decay of already produced ones (${\displaystyle \phi y}$, with defense persisting longer in the host when ${\displaystyle \phi }$ is decreased). The functions ${\displaystyle F}$ and ${\displaystyle G}$ describe the modulation of defense production, with ${\displaystyle F}$ modeling how production is activated upon pathogen detection, and ${\displaystyle G}$ describing negative regulation. Following Mayer et al., 1995, we wrote both ${\displaystyle F}$ and ${\displaystyle G}$ as non-linear functional responses and investigated how the shapes of these two functions influence the outcome of the infection.

The function ${\displaystyle F}$ is given by

We posited that ${\displaystyle F}$ must be an increasing function of pathogen load (${\displaystyle {\displaystyle u\ge v}}$ and ${\displaystyle {\displaystyle \beta >0}}$). In the absence of pathogens and with no negative regulation, the production of defense equals γ. The parameter γ, therefore, sets the maximum constitutive production of defense. Upon infection, ${\displaystyle F}$ increases with pathogen load (because we assume ${\displaystyle {\displaystyle \alpha >0}}$), and reaches its maximum feasible value when ${\displaystyle x}$ reaches the carrying capacity (i.e. when ${\displaystyle x=1}$ with then ${\displaystyle {\displaystyle F(1)=\gamma +\alpha /(1+\beta )}}$). Figure 2B illustrates that higher values of α, i.e., more sensitive defense activation, make this increase steeper. Increasing ${\displaystyle u}$ above ${\displaystyle v}$ makes ${\displaystyle F}$ sigmoid, which reproduces cases where the production of defense is barely activated unless the pathogen load becomes significant (as Figure 2B illustrates). Overall, the parametrization of ${\displaystyle F}$ allowed us to vary the constitutive level of defense expression, the speed at which defenses are activated upon infection, and the load at which pathogens are detected.

The function ${\displaystyle G}$ corresponds to the negative regulation of the production of defense, with

Here, we considered that the regulation of the production of defense must incorporate a negative effect of immune defense (${\displaystyle y}$) on defense production itself. This was motivated by the negative feedback loops that regulate immune responses in mammals and insects. In Drosophila, for example, the Diedel and WntD genes have such a function (Lamiable et al., 2016a; Lamiable et al., 2016b). Similarly, PGRP-LB, which has mammalian orthologs, is activated by the immune deficiency (Imd) pathway, one of the signaling cascades which regulates the expression of most antimicrobial peptides. But PGRP-LB also down-regulates this same Imd pathway, and thus provides a clear example of negative regulation of immune response (Kleino and Silverman, 2014).

Our model also aims to incorporate how damage accumulation modifies the regulation of defense production. This can conceptually be separated into two distinct phases. Early in the infection, damage-associated molecular patterns (DAMPs) have been proposed as universally used signals that trigger the activation of an immune response (Matzinger, 1994; Seong et al., 2021). These DAMPs can come from a wound, that would be systematically associated to the infection. We have assumed that in this case, DAMPs would trigger an instantaneous increase in both defense and damage, which we can model by increasing both ${\displaystyle y_0}$ and ${\displaystyle z_0}$. DAMPs could also originate from the infection itself, but they should then strongly correlate with bacterial load, as damage and load are tightly coupled during early infection. Therefore, we assumed that the positive regulation of defense production upon the detection of DAMPs that originate from the infection is incorporated in function ${\displaystyle F}$.

Later in the infection, the accumulation of damage caused by the infection could disrupt immune system homeostasis by altering the host’s health condition (Stoecklein et al., 2012). This may occur because the energy spent repairing damage is no longer available for defense production. Additionally, the tissues producing defenses may themselves be damaged by the infection, impairing the production of immune molecules or cells. In our model, these effects are represented in function ${\displaystyle G}$ by the negative effect of the damage (${\displaystyle z}$) on defense production (see Figure 2C), the importance of which is controlled by the parameters ${\displaystyle \psi }$ and ${\displaystyle l}$.

As should be clear from Equations 2 and 3, we purposely made ${\displaystyle F}$ and ${\displaystyle G}$ versatile so that our model can reproduce various types of host-pathogen interactions. This came at the cost of a large number of parameters. Figure 2 represents the model, lists its parameters and gives the parameter values we have used in our computations, unless otherwise specified. Parameters can be separated in two lists, with disease tolerance parameters that control damage production or repair on one side, and resistance parameters, which control the efficacy and level of immune defense on the other (see Figure 2). This clear-cut distinction, though, ignores potential indirect effects. For example, increasing the host capacity to repair damage (i.e. increasing ${\displaystyle \xi }$) makes it more tolerant to the infection. But damage hinders defense production, so that increasing ${\displaystyle \xi }$ indirectly allows the host to produce more defense, hence making it more resistant. Similarly, increased defense production causes additional damage and should, therefore, reduce the host’s tolerance to disease. One of our aims was, therefore, to describe how the combined direct and indirect effects of each parameter impact both resistance and tolerance to infections.

Finally, due to the complexity of functions ${\displaystyle F}$ and ${\displaystyle G}$, the equilibria of our model are intractable. We, therefore, relied on numerical methods to compute both the homeostatic state in the absence of pathogens ${\displaystyle ({\tilde{y}}_h,{\tilde{z}}_h)}$ and other equilibria ${\displaystyle (\tilde{x},\tilde{y},\tilde{z})}$ with pathogens present. This, as most of the numerical analysis detailed here, has been performed using the GSL library (Galassi et al., 2009) and functions provided in (R Developmental Core Team, 2018). We were nevertheless able to derive the stability conditions for some of the equilibria (see Appendix 1) and the relation between equilibrium load and parameters (see Appendix 2).

The entanglement between pathogen proliferation, defense, and damage production is too complex to predict the outcome of an infection from their separate effects. But numerically integrating Equations 1 allowed us to simulate infections, and thus to investigate how variations in parameters (e.g. the maximum level of constitutive production of defense, ${\displaystyle \gamma }$) or initial variable values (e.g. the initial level of damage, ${\displaystyle z_0}$) impact WHD and eventually determine the survival of the host.

By presenting a few chosen simulations, we first demonstrate that our model can reproduce situations where the host succeeds in controlling the infection (Figure 3A–C), either by clearing pathogens (${\displaystyle \gamma =5}$) or by maintaining pathogen load, and therefore damage, at a tolerable level (${\displaystyle \gamma =1.8}$). In the latter case, we can consider that the population of pathogens eventually reaches a Set Point Pathogen Load (SPPL, as defined in Figure 1) and will stay there for the whole host’s life.

Figure 3

Download asset Open asset

In the three simulations of Figure 3A–C, the level of defense was initially set at its homeostatic state, the immune state of a healthy host which is characterized by ${\displaystyle y={\tilde{y}}_h}$ and ${\displaystyle z={\tilde{z}}_h}$. Upon infection, defense increased and either returned to its initial level, when the infection is cleared, or stabilized to an intermediate level due to the remaining controlled pathogen population (Figure 3B, ${\displaystyle \gamma =1.8}$).

External challenges are expected to disrupt host homeostatic state prior to infection, which can lead to contrasted consequences. A host may be, for example, infected through a wound. This wound may cause damage which in turn could either trigger the immune system (Asri et al., 2019; Kenmoku et al., 2017), and potentially help the host to fight the infection, or hampers the expression of immune defense, and facilitates the infection (Stoecklein et al., 2012). We reproduced the later situation by setting ${\displaystyle z_0}$ above the homeostatic level of damage ${\displaystyle {\tilde{z}}_h}$. In our model, this means that defenses are initially hampered by the effects of damage, causing the immune response to occur with a delay (Figure 3B, gray curve with ${\displaystyle \gamma =5}$ and ${\displaystyle z_0={\tilde{z}}_h+10\mathrm{\%}{z}_d}$). This finding has been illustrated in Chambers et al., 2014, where the addition of sterile injury to Drosophila melanogaster prior to infection increased mortality to Providencia rettgeri compared to infections without wounding. This increase in mortality is suggested to be due to a decrease in resistance. In Figure 3B, the wound allowed pathogens to proliferate during early infection, but the host nevertheless managed to clear the infection. We demonstrated (see Appendix 1) that this is because the strong constitutive production of defense (${\displaystyle \gamma =5}$) in this simulation guaranties a high homeostatic level of defense. More generally, we proved that clearance is possible if and only if ${\displaystyle \gamma \ge {\gamma }_c}$, with

Infection resolution is thus possible only when a sufficient level of immune effectors is present before the infection starts. In our simulations, this level ${\displaystyle {\tilde{y}}_h}$ is set by the constitutive expression of defense ${\displaystyle \gamma }$, and characterizes the homeostatic state. As Equation 4 demonstrates, the level of constitutive expression required to cure the infection decreases when defenses are more efficient (${\displaystyle \delta }$ increases) or persist longer (${\displaystyle \phi }$ decreases), and increases with down-regulation of defense production.

The wound did not change the outcome of infection in hosts with strong constitutive immunity (gray curves in Figure 3A). It did in hosts with ${\displaystyle \gamma =1.8}$, which survived infection when not wounded (green curve in Figure 3A) but died from it when wounded prior to infection (${\displaystyle z_0={\tilde{z}}_h+0.3\mathrm{\%}{z}_d}$ in the green curves of Figure 3D). The initial immune handicap caused by the wound indeed facilitated infection, so that pathogen load reaches extreme values which eventually killed the host. Environmental challenges other than wounds might have a direct negative impact on defense production, which we could reproduce in our model by setting ${\displaystyle y_0}$ below ${\displaystyle {\tilde{y}}_h}$. This type of challenge could also facilitate infections by opportunistic microorganisms, with trajectories similar to the green curve of Figure 3D–F. As expected, hosts with weak constitutive defense can die from infection even when not wounded (${\displaystyle \gamma =0.5}$, red curves Figure 3D–F). In the two lethal infections presented in Figure 3D–F, the defense collapsed as damage rapidly amassed during the final phase of the disease (see Figure 3E), because ${\displaystyle \psi >0}$ makes damage hinder defense production. In both cases, damage eventually exceeded the maximum level the host can sustain (${\displaystyle z_d}$), and the host died from the infection (see Figure 3F).

Although identical in their final outcome, the two simulations of Figure 3D–F differ in the times the infection took to kill the host. When constitutive immunity is weak (red curve, ${\displaystyle \gamma =0.5}$), pathogens reached high loads before defense managed to slow down their proliferation. The hosts, therefore, succumbed rapidly. When constitutive immunity is stronger but the host is wounded (green curve, ${\displaystyle \gamma =1.8}$), the host went through a phase where pathogen proliferation was under temporary control, and infection transiently stabilized at a SPPL. This SPPL differs fundamentally from that of Figure 3A–C, because slowly accumulating damage ended up bringing defense below the level which permits efficient control of pathogen proliferation. This SPPL was, therefore, transient, while that of Figure 3A–C was stable.

The host depicted by the green curves of Figure 3A–C suffered from a chronic but benign infection; that of the green curves of Figure 3D–F died from a severe infection. Still these two hosts and their respective pathogens can be considered as genetically identical, as the two simulations were run with the very same parameter values. The only difference between them was that the host in Figure 3D–F had been wounded before being infected. This result demonstrates that our model can be bistable: for a fixed set of parameters, the outcome of an infection can depend on initial conditions.

We were not able to determine sufficient conditions which would guarantee bistability, but we demonstrated that bistability will occur if ${\displaystyle F}$ decreases for high pathogen loads (i.e. if ${\displaystyle u<v}$ and ${\displaystyle x}$ is large) or if ${\displaystyle G}$ decreases with damage (i.e. if ${\displaystyle \psi >0}$, see Appendix 1). We further demonstrated that this conclusion would hold if the damage repair mechanisms were regulated (see Appendix 1 for a complete analysis of bistability in modified versions of our model). The first condition for bistability, where ${\displaystyle F}$ decreases when ${\displaystyle x}$ gets large, implies that pathogens have a direct negative effect on defense production, as Ellner et al., 2021 assumed in their model. The second condition, where ${\displaystyle \psi }$ is non null, means that damage hinders the production of defense. We, therefore, propose that empirical evidence of bistability, such as in Figure 1, indicate that the infection lowers the immune defense either directly, as in Ellner et al., 2021, or indirectly, through resource diversion, as in van Leeuwen et al., 2019, or damage accumulation, as in our model.

The negative impact of infection on defense production is a necessary condition for bistability, but it is not sufficient on its own. Additional, more specific conditions are required, which cannot be mathematically determined. In the following, we focused on two parameters which describes the constitutive and the inducible parts of defense production (namely γ, the maximum constitutive defense production and α, the activation of defense production). In Figure 4A, where the two parameters were varied, two lines delimit a region within which the infection dynamics is bistable. When γ is fixed to 2 and only α is varied (which corresponds to the vertical dotted line in Figure 4A) the equilibrium load ${\displaystyle \tilde{x}}$ is very high for low values of α and drops to low values when α is high (Figure 4B). For intermediate values of α, the system is bistable: both high load and low load equilibria are possible. The same sort of pattern occurs when α is fixed to 10 and γ only is varied (the horizontal dotted line in Figure 4A): the equilibrium load is high when γ is low, the infection is always cleared when ${\displaystyle \gamma }$ is high, and the system becomes bistable when γ has intermediate values. Inside this bistable region, the high load infection is always possible; the second possible equilibrium is either low load infection, when ${\displaystyle \gamma <{\gamma }_c}$, or clearance otherwise (Figure 4B).

Figure 4

Download asset Open asset

In bistable situations, the equilibrium load is either high or low, with no possible intermediate situation. This happens because in such systems the two possible stable equilibria are separated by a third, unstable, one that the system cannot reach (the dashed curves in Figure 4B and C). Which stable equilibrium will be reached then depends on initial conditions. This is represented by the colored areas of Figure 4A, with the inoculum ${\displaystyle x_0}$ fixed to 10⁻⁶ and ${\displaystyle y_0}$ and ${\displaystyle z_0}$ set to the homeostatic state. When only the high load equilibrium is possible, the host always dies from the infection because the equilibrium level of damage ${\displaystyle \tilde{z}}$ exceeds ${\displaystyle z_d}$. When only the low load equilibrium is possible, the host survives, and the infection becomes chronic. In the bistable area, the infection also becomes chronic when clearance is not possible (${\displaystyle \gamma <{\gamma }_c}$), with load stabilizing at a SPPL. But two situations must be distinguished here: first, when the immunity is strong enough to control pathogen’s proliferation (the green area in 4 A), the infection reaches the low load equilibrium. The host will then survive the infection and load will stabilize at the SPPL for its whole life. Second, when the host immunity is weaker (the yellow area in Figure 4A) defense controls proliferation but for a limited period of time, as in the green curve of Figure 3D. The SPPL is then transient: the infection will eventually break loose, as accumulating damage causes the collapse of defense, ultimately killing the host.

Demonstrating experimentally that a SPPL is transient might prove difficult, first because it is not always feasible to monitor pathogen load until the host dies and, second, because even if it was, it may still be challenging to demonstrate that the death of an individual host was actually caused by the infection. Our model offers a solution to this problem. Mathematically, a transient SPPL indeed occurs when the three variables progress slowly, because they have approached an equilibrium point, where all derivatives ${\displaystyle dx/dt}$, ${\displaystyle dy/dt}$ and ${\displaystyle dz/dt}$ are close to zero. However, this equilibrium must be unstable, as the trajectories eventually moves away from it and head towards a high load infection. This indicates, first, that the SPPL cannot be transient unless an unstable equilibrium point exists with ${\displaystyle \tilde{x}>0}$, which in turn requires that the system is bistable (this unstable equilibrium corresponds to the dashed curves in Figure 4B and C). Second, it also proves that a transient SPPL is not a true equilibrium situation. It is in fact close to the notion of a quasi-static state in thermodynamics. Therefore, the transient SPPL should vary with initial conditions. We have shown that the SPPL does indeed increase with the initial dose ${\displaystyle x_0}$ when it is transient, while it is independent from initial conditions when stable (see Figure 5). These results are in line with empirical results, as the SPPL depends on initial inoculum when Drosophila are infected by pathogenic bacteria (Chambers et al., 2019; Duneau et al., 2017a; Acuña Hidalgo et al., 2021), but not when infected by non-pathogenic bacteria such as E. coli (Ramirez-Corona et al., 2021). We, therefore, propose that testing the link between the SPPL and the initial dose is a way to experimentally demonstrate that the control over infection is transient.

Figure 5

Download asset Open asset

Bistability, as we have already discussed, is possible in our model when accumulating damage negatively impacts defense production. As expected, reducing this impact by lowering the value of ${\displaystyle \psi }$ reduces the range of parameters which permits bistability. With ${\displaystyle \psi =2}$ in Figure 4A, 94% of the parameter area where the high load equilibrium is possible is bistable; in Figure 4D, where ${\displaystyle \psi }$ is reduced to 0.5, only 82% of this same area is bistable. Finally, the distinction between high and low load equilibria, although convenient, is not always possible. In Figure 4E, with ${\displaystyle u=0.75}$, the bistable area forms a closed triangular shape. In this case, when low constitutive defense production (i.e. low γ) does not allow for bistability, a continuum exists between high-load equilibria (low values of α) and low-load equilibria (high values of ${\displaystyle \alpha }$), without a bistable region in between.

The SPPL and the PLUD can both be estimated from experimental infections, as earlier work by Duneau et al., 2017a has shown. They have been proposed to reflect different aspects of the host immune response: the SPPL, stable or transient, being a pathogen load stabilized by the host immunity, may be taken as a proxy for resistance; the PLUD, being the maximum pathogen load the host can endure, may rather quantify the host’s tolerance to infection. Measurements of host survival to infection are also commonly used as proxies of resistance or tolerance (e.g. Gupta and Vale, 2017; Louie et al., 2016). Here, we have used our model to question the way these quantities are used as surrogate measurements for tolerance or resistance. For this purpose, we computed them from simulated infections and surveyed how they relate to immunity parameters (see Appendix 3,4 and 5 for a complete analysis).

As expected, the time it takes for a pathogen to kill its host decreases if inoculum size is increased (Figure 5A) or if the host is wounded prior to infection (i.e. when ${\displaystyle z_0}$ is increased by 1% of ${\displaystyle z_d}$ above the homeostatic state). Lethal Time (LT) thus strongly depends on the conditions in which the infection has been initiated. Our simulations further show that any genetic change in the host or the pathogen that would reduce the host resistance (red curves in Figure 5A) or tolerance (blue curves in Figure 5A) would accelerate death. We also showed that contrary to what some authors have proposed (e.g. Gupta and Vale, 2017; Louie et al., 2016), a reduction in tolerance does not necessarily modify how LT relates to dose; a reduction in resistance, conversely, makes LT less dependent on dose (compare red and blue curves in Figure 5A, see Appendix 5 for a complete analysis on this point). Hence, the comparison of lethal time in response to different doses is more likely to characterize differences in resistance than in disease tolerance.

The SPPL, when it is stable, does not depend on initial conditions (dashed lines in Figure 5B). This is because a stable SPPL is an equilibrium towards which the infection will converge, no matter where it starts from, as long as the host survives the infection. The transient SPPL, conversely, is not an equilibrium. We showed that it does increase with dose and is higher when the host is wounded prior to infection. Our simulations finally indicate that both transient and stable SPPL increase when either tolerance or resistance are impeded (see Appendix 1 and 3 for a complete analysis). In fact, if by definition stable SPPL does not correlate with the time the host will die, variations in the transient SPPL almost exactly mirrors that of the LT. This is because the control of proliferation does not last long when the transient SPPL is high: infections which maintain high pathogen loads rapidly progress towards host death. A first practical consequence of this observation is that measurements of the SPPL and of the LT bear the same information on the host immunity. Therefore, it is useless to measure both, unless one wants to predict transmission rate, which should quantitatively depend on both bacterial load and infection duration. A second practical consequence is that the SPPL, when transient, can be used as a predictor of the infected host’s lifespan, as already demonstrated in HIV infections (Fraser et al., 2007; Mellors et al., 1996).

The two hosts of Figure 3D–F die at very different times but at almost indistinguishable PLUD. This is confirmed by Figure 5A, where the PLUD does not vary with changes in inoculum size that do otherwise yield threefold variations in LT (compare Figure 5C to Figure 5A). More generally, the PLUD is only weakly influenced by the conditions in which the infection has been initiated, albeit we predict that it should slightly increase when the host is wounded. Our model, therefore, reproduces an important property of the PLUD which Duneau et al., 2017a have documented (see Figure 1): for a given pair of host and bacterial pathogens, the PLUD in D. melanogaster is almost constant and does not correlate with the time to succumb from the infection. In our model, this comes from the fact that the pathogen load ${\displaystyle x}$ evolves much faster than other variables at the start of the infection (as detailed in Appendix 4). Pathogens thus rapidly reach the highest possible load allowed by the current amount of defense (${\displaystyle x_t\lesssim 1-\delta y_t}$). The dynamics of ${\displaystyle x}$ is then ‘enslaved’ by that of ${\displaystyle y}$ and ${\displaystyle z}$ (Van Kampen, 1985) and, as the disease enters its final stage, the load increases slowly while ${\displaystyle y}$ gradually diminishes under the effect of accumulating damage. Because of this very particular dynamics, the way ${\displaystyle x_t}$ relates to ${\displaystyle z_t}$ will be very similar for all infections that enter their final phase, which in the end renders the PLUD almost independent from the conditions that prevailed at the onset of the infection.

Another particularity of the PLUD, is that it increases when resistance is diminished (red curves in Figure 5C) but decreases when tolerance is lowered (blue curves in Figure 5C). The PLUD is, therefore, a composite measurement which reacts to both resistance and tolerance traits, just like the SPPL and the LT, but has the unique feature that it responds differently to these two types of variation (see Appendix 4 and the following section for a more complete elasticity analysis).

We have seen before that variations in resistance and variations in disease tolerance have distinctive impacts on the PLUD. Still, a host lineage with higher than average PLUD could either be less resistant or more tolerant than other lineages. The PLUD is, therefore, difficult to interpret by itself, but we shall demonstrate that its relationship with mortality measurements are informative. To investigate this, we performed elasticity analyses of the PLUD and the Hazard Ratio, a measure of death risk relative to a chosen reference host (HR, see Appendix 6 for definitions and analyses). As our model allows to determine the time to death of a given infection, the HR can be obtained from a Cox proportional hazard model (Cox, 1972) adjusted on simulated survival curves.

As expected, we found that reducing the risk to die from the infection can be achieved through a tighter control of pathogen proliferation or, alternatively, by mitigating the damage the infection causes (Figure 6A and B). For example, HR decreases both when the defense is made more efficient (${\displaystyle \delta }$ increases) and when damage repairing is accelerated (${\displaystyle \xi }$ increases).

Figure 6

Download asset Open asset

Variations in PLUD are more complex, but Figure 6A suggests that reducing damage (e.g. by increasing damage repair, ${\displaystyle \xi }$) or making the host more damage-tolerant (by increasing the maximum damage level the host can sustain, ${\displaystyle z_d}$) should increase the PLUD. As a result, variations in parameters which have a direct impact on the dynamics of damage produce a negative correlation between PLUD and HR. Figure 6B further suggests that anything that reduces pathogen proliferation, like increasing the production or efficacy of defense (α or ${\displaystyle \delta }$), comes with a decreased PLUD. Variations in parameters that control the defense, therefore, result in a positive correlation between the PLUD and HR. In summary, Figure 6 indicates that faster death (i.e. increased HR) would indicate lower tolerance when associated with lower PLUD, when it would rather suggest lower resistance in case of higher PLUD.

The case of ${\displaystyle \eta }$, which quantifies the importance of defense side effects, demonstrates that although appealing, this method has its limits. This parameter clearly relates to damage production and we demonstrate that the lower bound of the PLUD expectedly decreases when ${\displaystyle \eta }$ increases (see Equation (S4-1) in Appendix 4). But when immunity is weak, this effect can reverse, with increasing ${\displaystyle \eta }$ impairing control and resulting in higher PLUD (as we demonstrate in our elasticity analysis, see Appendix 4—figure 2 when both α and γ are low). In hosts with very weak immunity, therefore, variations in the detrimental side-effects of defense (${\displaystyle \eta }$) can create a positive correlation between PLUD and HR (as in Figure 6A). In addition, Figure 6A shows that different forms of variation in tolerance (e.g. variations in pathogenicity ${\displaystyle \omega }$, in damage repair ${\displaystyle \xi }$ or in ${\displaystyle z_d}$ the maximum level of damage the host can sustain) will produce the same signal. Figure 6B demonstrates a similar situation for variations in resistance. Therefore, experimental measurements of PLUD and HR can be used to distinguish variations in tolerance from variations in resistance, broadly speaking, but will likely fail to identify the precise mechanisms involved.

The practical use of the above proposed method will largely depend on which specific parameters differ the most among the hosts and pathogens to compare. We cannot predict much on this matter, as everything will obviously depend on the genetic variation of these traits. A crude way to test the utility of the method, though, is to let all parameters vary and see which effect predominates. In Figure 6C, we randomly drew all 14 parameters from independent Gaussian laws with average set at the parameter values given in Figure 2, and a variance chosen so that 95% of the random values differ by less than 0.5% from this average. If we consider that the parameter set given in Figure 2 represents a wild-type host, each random parameter set can be considered as being a mutant strain. We randomly sampled 2000 such mutants. We then simulated 100 survival curves for each of them, which we compared to that of the wild type by using a Cox proportional hazard model to estimate HR. We used the same simulations to compare the PLUD of the mutants to that of the wild-type, as if load was estimated by plating. For that purpose, we randomly drew numbers of Colony Forming Units (CFU) from Poisson distributions with average given by the predicted WHD, which we then compared between the mutant and the wild-type hosts using a Poisson glm. We have drawn one number of CFU per host from a Poisson distribution with average 1000 times the simulated PLUD, so that we have similar statistical powers when comparing PLUDs and HRs. In Figure 6C, the red cross lying on the origin is the wild-type; mutants which significantly deviate from it (closed circle) are separated in four groups, according to the deviation signs. For each group, we compared the distribution of mutant traits to that of the wild-type. We found that mutants with significantly higher HR and PLUD have almost systematically a stronger down-regulation of defense production (with clearly higher value of ${\displaystyle l}$) and a slightly lower efficiency of defense (${\displaystyle \delta }$). Hence, such mutants die faster than the wild-type because they are less resistant. Mutants with higher HR but lower PLUD have most often a lower tolerance to damage (lower ${\displaystyle z_d}$) and are more susceptible to the pathogen virulence (greater ${\displaystyle \omega }$). These mutants, therefore, die faster than the wild-type strain because they are less disease tolerant.

Our model predicts that when damage hinders defense production, hosts that are wounded prior to infection should have higher PLUD (see Figure 5C). We tested this prediction by measuring the PLUD of the bacterial pathogen Providencia rettgeri injected in D. melanogaster that were either injured in the thorax before being injected in the abdomen or not injured (see Figure 7 and Appendix 7 for methodological details). Wound could in principle increase defense, e.g. as hemocyte innate immune training triggered by DAMPs (Chakrabarti and Visweswariah, 2020), but hemocytes have no detectable effects on infections caused by Providencia rettgeri (Duneau et al., 2017a). Therefore, we do not anticipate any positive impacts of the wound. Instead, we expect that the wound decreases resistance, as indicated by Chambers et al., 2014. We found that indeed, the wound significantly reduced survival (Cox proportional hazard model with random experimental blocks: X² = 1397.44, degrees of freedom (df) = 1, p - value < 2.2e^-16). PLUD estimates were found to be highly variable, with some individuals dying at load less than 10⁴. These individuals may have died from reasons other than the P. rettgeri infection; we, therefore, conducted two analyses, with these individuals either included or removed (outliers being identified by a Rosner test, see Appendix 7 for details). We found in both analyses that wounded hosts die at a PLUD significantly higher than non-wounded (linear mixed model with random experimental blocks: X² = 11.29, df = 1, p - value = 0.008 with outliers removed; X² = 11.04, df = 1, p - value = 0.0009 with outliers retained) with no significant difference among genotypes (p - value = 0.73 and p - value = 0.59) with outliers removed or not removed, respectively although the effect is much clearer on line RAL-630 than in the three other lines (see Figure 7). Overall the wound does increase both hazard ratio and PLUD, as our model predicts. This result can be taken as an evidence that inflicting damage to D. melanogaster hinders the production of defense, and thus indirectly reduces resistance.

Figure 7

Download asset Open asset

Our model also predicts that mutations that decrease resistance should in most cases increase the PLUD. This seems in contradiction with Duneau et al., 2017a who did not find any effect of immunosuppression on the PLUD in P. rettgeri infections. The original analysis in this paper used a linear mixed model where the effects of Imd, Melanization and Toll pathways were assessed together. We reanalized the dataset and examined each pathway separately. We observed that the PLUD exhibited an increase when the Imd pathway was suppressed (Kruskal-Wallis: X² = 7.87, df = 1, p-value = 0.005) or when melanization was inhibited (Kruskal-Wallis: X² = 15.56, df = 1, p-value = 7.9e^-5). However, there was no significant change in PLUD when the Toll pathway was suppressed (Kruskal-Wallis: X² = 0.03, df = 1, p-value = 0.86).

A limitation of the previous study is that the wildtypes used were not ideal genetic controls, given that the mutants did not undergo backcrossing. In order to strengthen the test of our theoretical prediction, we quantified the PLUD in D. melanogaster mutants with some immune effectors (AMPs) either deleted or silenced, and compared it to PLUD measured in appropriate controls (see Appendix 7).

In a first experiment, we used two isogenic lines sharing the same genetic background but from which genetic deletions have induced a loss-of-function either in Defensin (line A) or in Drosocin, Attacins, and both Diptericins (line B). We also used a third mutant line (line AB) which combines the loss-of-function of both lines A and B. Line A lacks Defensin, an antimicrobial peptide strongly upregulated upon P. rettgeri infection (Duneau et al., 2017b; Troha et al., 2018) but with moderate effects on survival (Hanson et al., 2019). Lines B and AB, which both lack Diptericins, are conversely expected to be highly susceptible to P. rettgeri (Hanson et al., 2019; Unckless et al., 2016). We have used as a control a fourth mutant line which lacks Bomanins (line Bom^Δ55C). Bom^Δ55C is a good reference as it is mutated in the same background than the other genotypes, but the AMPs it lacks are mostly involved into fighting Gram-positive bacterial infections (Clemmons et al., 2015). In addition, previous experiments have shown that Bom^Δ55C mutants have a PLUD comparable to controls when infected with the gram positive bacteria Enterococcus faecalis (Lin et al., 2020).

We found that, compared to line Bom^Δ55C, lines A, B, and AB have higher death rates when infected by P. rettgeri (Figure 8A; Cox proportional hazard model with random experimental blocks: X² = 138.12, df = 3, p-value = 9.6e^-30) and higher PLUD (Figure 8B; linear mixed model with random experimental blocks: X² = 17.01, df = 3 , p-value = 0.0007). We challenged this first result by comparing another mutant that, contrary to line B, lacks only Diptericins (Dpt^SK1) (Hanson et al., 2019). We compared this mutant to Drs^R1, a mutant which shares the same genetic background (Hanson et al., 2019) but has Drosomycin deleted, Drosomycin being an antifungi AMP that is inactive against bacteria (Fehlbaum et al., 1994). We found again that death is faster and PLUD higher when Diptericins are deleted (Figure 8D–F; linear mixed model with random experimental blocks: X² = 10.11, df = 1, p-value = 0.001). But if Bom^Δ55C and Drosomycin are clearly less important for resistance against P. rettgeri than Diptericin, studies start to suggest that they may play a role in tolerance to fungi toxin or are involved in other functions which could influence tolerance (Araki et al., 2019; Xu et al., 2023). Therefore, we also challenged our hypothesis by injecting P. rettgeri in flies that constitutively express, thanks to an ubiquitous driver (Actin5C-Gal4), a RNAi which silences specifically Diptericin A. This time we used as control line the advised match genetic background for the attP40 genetic background of the TRiP genetic RNAi panel. We found that flies with silenced Diptericin A die faster and at a higher PLUD than controls (Figure 8G–I; linear mixed model with random experimental blocks: X² = 8.42, df = 1, p-value = 0.0037). Overall, based on previous results and on those three different approaches, we concluded that a reduction in resistance caused by a lack of effectors provokes an increase in both HR and PLUD. Finally, it should be noticed that, as our model predicts, the increase in PLUD is remarkably constant throughout experiments (approx. +0.2 in Figure 8C, F and I) even though the HR varies from approx. +0.3 to +1.5.

Figure 8

Download asset Open asset

No gene is unambiguously identified as involved in disease tolerance in D. melanogaster, probably in part because damage are uneasy to quantify experimentally. CrebA and Bombardier are candidate genes (Troha et al., 2018; Lin et al., 2020), but they have been determined as such using the PLUD as a tolerance proxy, which would make our reasoning circular. Other genes have been proposed, but their roles were mostly identified using approaches that our model suggests are limited for establishing a proxy of disease tolerance (Figure 5A, see Appendix 5). We thus decided to test our method on a gene which function is well understood and should contribute to disease tolerance.

Reactive oxygen species (ROS) have been shown to be critical agents in oxygen toxicity, disrupting the structural and functional integrity of cells. Catalase is among several enzymes involved in scavenging oxygen free radicals and protecting those cells. The activity of Catalases into reducing reactive oxygen species has been shown in Drosophila infections (Ha et al., 2005) and mosquito systemic infections (DeJong et al., 2007), making it a strong candidate for a tolerance gene.

We found that silencing Catalase ubiquitously with Actin-Gal4 driving an RNAi does increase flies susceptibility to systemic bacterial infection when comparing to a match genetic background as control (i.e. attP2 control line, the advised RNAi control for the attP2 genetic background of the TRiP genetic RNAi panel, see Figure 9A). The log(HR) is comparable in strength to that of silencing Diptericin A. We further observed that silencing Catalase does decrease the PLUD (see Figure 9B; linear mixed model with random experimental blocks and a correction for heteroscedasticity: X² = 23.33, df = 1, p-value = 1.37e^-6). Again, although the level of PLUD was different, likely due to the difference in genetic background, the relative difference (i.e. log(PLUD ratio)) was comparable in strength to that of silencing Diptericin A, but of opposite sign. Therefore, based on its known functions and the fact that our experimental results follow the predictions of our theoretical model, we propose that Catalase is part of the arsenal that enable a fly to tolerate bacterial infections.

Figure 9

Download asset Open asset

When a host survives the infection but does not successfully clear the pathogen, the disease enters a phase of chronicity (Virgin et al., 2009) and load stabilizes at a SPPL. Such chronic phases characterize a large range of infectious diseases, and we have shown that this is determined by the host immune defense. Notably, strong constitutive defense prevents chronicity because it permits the host to clear pathogens; weak constitutive defense combined to slow activation upon pathogen detection also prevents chronicity because hosts then die from the infection.

A more surprising prediction of our model is that chronicity can be transient, the SPPL being then unstable. This happens when the damage which accumulates during the infection eventually obstructs the production of immune defense. Pathogen proliferation is then unleashed which leads to host death. Illustrations of such infections could be, in humans, Mycobacterium tuberculosis causing tuberculosis infections or HIV causing AIDS; in Drosophila, this could be Providencia rettgeri infections. When the SPPL is stable, by contrast, the host never clears pathogens but it will not die either from the infection. In human, this could be Porphyromonas gingivilis causing gum infections; in Drosophila, this could be Escherichia coli infections. Our model clearly highlighted those two types of chronic infections and suggests that their SPPL cannot be interpreted the same way.

In particular, we found that unstable SPPLs depend on how the infection was initiated. For example, we predict that they should be higher when the host is injured prior to contamination and should correlate positively with the inoculum size (see Figure 5). These predictions agree with previous study of bacterial infection of D. melanogaster (Acuña Hidalgo et al., 2021; Chambers et al., 2014; Chambers et al., 2019; Duneau et al., 2017a). Our model predicts that stable SPPLs are, conversely, independent from how the infection is initiated. In particular, they are insensitive to inoculum size, in accordance with benign Escherichia coli infections of D. melanogaster (Ramirez-Corona et al., 2021). The stable SPPL, therefore, reflects stable genetically determined characteristics of both the host and the pathogen, but it is affected by traits related to both tolerance and resistance. For example, low stable SPPL could indicate either that the host has efficient immune defense, or that it tolerates well infections.

This is not to say that no information can be gained from the study of SPPL, be it stable or unstable. In probably all chronic infections, everything being equal, the pathogen load during the chronic phase is an appropriate tool to compare the capacity to be transmitted to another host: the higher the SPPL, the most the host sheds pathogens (Fraser et al., 2007; Matthews et al., 2006). The unstable SPPL is also an important concept because it is a good predictor of the duration of chronicity: we found indeed that the higher the SPPL, the earlier the host dies from the infection. Monitoring survival to infection and transient SPPL is, therefore, providing redundant information. This is particularly well established in the specific case of HIV infections, where high SPPL is associated to short asymptomatic phases and rapid progression to the final stage of the disease (Fraser et al., 2007). The SPPL in HIV infections has also raised considerable attention because of two striking characteristics which are not fully understood. First, it is highly variable among patients; second, it is heritable, which means that the SPPL of a newly infected patient positively correlates to that of the person it has been infected by. The most common explanation of this correlation is that part of the variation in SPPL is caused by genetic variation in viruses. But if variance in SPPL has a genetic origin, the SPPL should evolve during the course of the infection and load should, therefore, increase over time.

The heritable nature of SPPL, therefore, contradicts the observation that load is constant throughout the asymptomatic phase of the disease (although some possible solutions to this paradox have been proposed Bonhoeffer et al., 2015; Hool et al., 2014). Our model suggests that part of the large variance in SPPL could be non-genetic yet heritable. If transmission events occur during chronicity, the inoculum size should indeed increase with the SPPL of the donor host (as suggested by Fraser et al., 2007). As we have demonstrated that the transient SPPL does increase with inoculum size, the SPPL of a newly contaminated host should reflect that of the donor host. The fact that the transient SPPL increases with inoculum size could, therefore, suffice to create a form of non-genetic heritability in SPPL.

Finally, we propose that quantifying the SPPL for various doses in a chronic disease would allow, if the SPPL increases with dose, to demonstrate that it is unstable. This could be used as an experimental method to prove that chronicity is transient. In experimental systems which our model reproduces, it would in addition demonstrate that the system is bistable (a necessary condition for the existence of unstable SPPL) which in turn requires that the infection impairs immunity by reducing the production of defense, either directly as in Ellner et al., 2021 or indirectly through the accumulation of damage (Stoecklein et al., 2012; van Leeuwen et al., 2019; DeJong et al., 2007), or by making defense less efficient (e.g. Zhang, 2016). Hence, showing that the SPPL is transient informs on the fact that the pathogens, or the damage they cause, obstruct the immune response. The positive link between dose and SPPL has been found in some experiments (Shultz et al., 2016; Chambers et al., 2019); additional work is required to demonstrate that the explanation our model provides applies to these experiments.

We have so far assumed that pathogens do not evolve during the course of the infection. We did not consider this possibility, first because it is beyond the scope of our study, and second because Duneau et al., 2017a demonstrated that within-host evolution did not explain bifurcation in their experiments. But in situations where, for example, frequent mutations arise that make pathogens resistant to host immunity, a phenomenon resembling bistability could occur: hosts where resistance has evolved would be killed by the infection while others would control it and survive. If this happens, bistability would have nothing to do with damage hindering defense. But this does not compromise our predictions because the SPPL would then probably not vary with dose (as demonstrated by Ellner et al., 2021 in the modified version of their model where bacteria are protected from AMPs).

The Pathogen Load Upon Death is the maximum pathogen load a host can sustain before it dies from the infection. It has been logically used as a measure of disease tolerance (Duneau et al., 2017a; Huang et al., 2020; Lin et al., 2020; Troha et al., 2018) and of pathogenicity (Faucher et al., 2021). As expected, our model confirms that the PLUD strongly depends on tolerance to damage (${\displaystyle z_d}$, see Appendix 4). However, if the PLUD was an unambiguous quantification of tolerance, as proposed by Duneau et al., 2017a, it should not depend on any parameter other than ${\displaystyle z_d}$, or at least not on parameters that determine the control of pathogen proliferation. We conversely predict that hosts with an impeded immune response are expected to have a higher PLUD, which we confirmed with an experimental design more powerful than that of Duneau et al., 2017a (see Figure 8).

In our model, the PLUD increases when less defense are produced (e.g. when α or γ decrease, see Appendix 4—figure 2) or when defense are made less efficient (i.e. when ${\displaystyle \delta }$ is decreased, see Figure 6). This is because the PLUD is a load which, as such, reflects the ability of pathogens to proliferate inside the host. Hence, like the transient SPPL and probably like any quantity derived from load estimations, the PLUD is subject to mixed influences of immune control of proliferation and of damage mitigation.

The PLUD, unlike the transient SPPL this time, is independent from inoculum size: flies injected with large inoculum die faster but at the same load than flies infected with small inoculum (Duneau et al., 2017a). Our model reproduces this observation, and we have demonstrated mathematically that this happens because pathogen load is a ‘fast variable’ (Van Kampen, 1985, see Appendix 4). In practice, once the accumulation of damage due to the infection is starting to affect the immune system in a way that the defense are impeded, the system is entering a run-away process: pathogens proliferate faster, inflicting new damage to the host, which decrease the defense further. Infections that enter this dynamics tend to all follow the same typical trajectory, and therefore kill the host at the same pathogen load. We found that this occurs in most parameter combinations (see Appendix 4), except when hosts die at very low damage (i.e. when ${\displaystyle z_d}$ is low). In that particular case, death occurs so rapidly that infections have no time to converge to their typical final trajectory and the PLUD depends on inoculum size. This could represents situations where the host is already extremely weak or about to die at the time of the infection. However, in most relevant cases, the PLUD being independent of how the infection has started, it reflects the genetically fixed traits of both hosts and pathogens that influence disease tolerance and resistance.

It has been proposed that disease tolerance could be quantified as the slope of a linear relation of time to death to a load measured at a fixed point during the infection (e.g. Louie et al., 2016; McCarville and Ayres, 2018) or with the linear or logistic relationship between initial inoculum size and time to death (e.g. Gupta and Vale, 2017). However, our analysis strongly suggests that these approaches mostly measure resistance (see Figure 5). This is probably due to the fact that any pathogen load reflects resistance mechanisms. Instead, we propose to use the relation between the HR and the PLUD ratio to determine whether differences in susceptibilities are mainly due to deficiencies in tolerance or resistance mechanisms. Our method is similar but differs in two key aspects. First, we propose using a load that reflects a specific moment of the infection rather than a load measured at an arbitrary, fixed time. This approach ensures that loads can be compared among hosts and among pathogens, as they will always have the same biological interpretation. Second, we used measures of Hazard Ratio obtained for a fixed inoculum size, as changing the inoculum size has little impact on the relation between the PLUD and the Hazard Ratio, unless very high inoculum are used, which can lead to artifacts by artificially bypassing the immune response.

Observing that a particular genotype of host has a lower than average PLUD cannot indicate whether this genotype is less tolerant or more resistant than others to the infection. However, our model shows that if this low PLUD genotype also has a higher than average risk to die from the infection (or a shorter time to death) its susceptibility is probably due to a default in tolerance to infection. If conversely a genotype has a higher than average risk to die with a higher PLUD, we would conclude that its susceptibility is mostly due to a default in resistance (see Figure 6).

We confirmed this prediction in the case of P. rettgeri infecting D. melanogaster: host genotypes with impeded immune effectors have both a higher Hazard Ratio and a higher PLUD than host genotypes with a more efficient immune response. Using Catalase, an enzyme involved in protecting cells from damage by ROS, as a candidate gene for tolerance, we also confirmed that faster death, which we assumed to indicate lower tolerance, was associated with a lower PLUD. Although the link between a lower pathogen load sustained before death and a default in tolerance seems logical, this probably needs further confirmation with other candidate genes for which we will start to understand the mechanisms by which they contribute to tolerance. As we already mentioned, potential tolerance genes in D. melanogaster have been identified using the PLUD as a tolerance proxy (CrebA and Bombardier Lin et al., 2020; Troha et al., 2018), or methods of which we discussed the limits, (grainyhead, debris buster, and G9a Howick and Lazzaro, 2017; Merkling et al., 2015). In fact, the role of G9a in controlling bacterial load is still debated (Gupta and Vale, 2017). We, therefore, here provide a clear indication that our method allows to detect defaults in tolerance, but complete validation will require that other candidate genes are tested.

One of the great benefit of theoretical studies is that they allow to explore notions which arose from empirical work. Concepts which seem well and clearly defined may often be underlain by complex and intermingled mechanisms. Our approach has helped to clarify the notions of disease tolerance and resistance, in particular by showing how experimental measurements can be misleading when interpreted as unambiguous measurements of one or the other mechanism. Many other mechanisms which determine important aspects of infection pose the same sort of difficulty. For example, the co-regulation of damage repair and immune response is now well accepted (Martins et al., 2019; Pucillo and Vitale, 2020) but remains difficult to study. Similarly, we do not know how the presence of different types of damage could affect disease outcome. We advocate that pursuing this combination of theoretical and empirical approaches is needed, first to predict how such intermingled mechanisms impact the infection outcomes, and second to design appropriate proxies and experimental methods to quantify their impact.

Let ${\displaystyle \tilde{x}}$, ${\displaystyle \tilde{y}}$, and ${\displaystyle \tilde{z}}$ be the state variables at equilibrium. Let ${\displaystyle \overline{y}(x)}$ and ${\displaystyle \overline{z}(x)}$ be the ${\displaystyle y}$ and ${\displaystyle z}$ values such that ${\displaystyle \frac{dx}{dt}(x,\overline{y}(x),\overline{z}(x))=\frac{dz}{dt}(x,\overline{y}(x),\overline{z}(x))=0}$. From Equation 1, it is easily shown that

with ${\displaystyle \tilde{y}=\overline{y}(\tilde{x})}$ and ${\displaystyle \tilde{z}=\overline{z}(\tilde{x})}$. The equilibrium load ${\displaystyle \tilde{x}}$ can be found by solving ${\displaystyle dy/dt(x,\overline{y}(x),\overline{z}(x))=0}$, which can be written as

The system is, therefore, bistable only if ${\displaystyle H}$ has at least two roots in [0, 1], which implies that ${\displaystyle H^{\mathrm{\prime }}}$ has at least one root in the same interval. We obtain

with ${\displaystyle \dot{\overline{y}}=d\overline{y}/dx=-1/\delta }$ and ${\displaystyle \dot{\overline{z}}=d\overline{z}/dx=(\omega -\eta /\delta )/\xi =\theta /\xi }$ is positive under the assumption that ${\displaystyle u\ge v}$, and all derivatives of ${\displaystyle G}$ are negative. If ${\displaystyle \psi =0}$, ${\displaystyle \mathrm{\partial }G/\mathrm{\partial }z=0}$ and all remaining terms in Equation S1-6 are positive. ${\displaystyle H^{\mathrm{\prime }}}$, therefore, never cancels and ${\displaystyle H}$ must have a single root: the system cannot be bistable. If, conversely, ${\displaystyle \psi >0}$ and if ${\displaystyle \theta =\omega -\eta /\delta >0}$, then ${\displaystyle \dot{\overline{z}}\mathrm{\partial }G/\mathrm{\partial }z<0}$ and ${\displaystyle H}$ can have more than one root. ${\displaystyle \psi >0}$ and ${\displaystyle \theta >0}$ are, therefore, necessary conditions for the system to be bistable.

We shall now investigate the robustness of this conclusion when three assumptions of our model are relaxed. We have assumed so far that

• (i) Damage repair occurs at a constant rate ${\displaystyle \xi }$. Relaxing this assumption would make ${\displaystyle \xi }$ (or even ${\displaystyle \omega }$ and ${\displaystyle \eta }$) depend on ${\displaystyle x}$. The derivative ${\displaystyle \dot{\overline{z}}(x)}$ would then be modified,

• (ii) Damage accumulation can only decrease defense production. But it is known that host cells attacked by pathogens can emit molecules which trigger immune defense. When this occurs, defense production may increase with damage, at least while damage stays moderate. This situation can be reproduced in our model by letting ${\displaystyle F}$ increase with ${\displaystyle z}$,

• (iii) Defense are not consumed when they kill pathogens. Many other models (e.g. Ellner et al., 2021; Souto-Maior et al., 2018) rather include a term which reduces the amount of defense in proportion to their anti-pathogen activity (i.e. they include a negative term proportional to ${\displaystyle \delta xy}$).

If only the first assumption is relaxed, the derivative ${\displaystyle \dot{\overline{z}}(x)}$ is changed, as we said, but it will still cancel out in Equation S1-6 when ${\displaystyle \psi =0}$. Modulation of damage repair alone would, therefore, not invalidate our conclusion that bistability is impossible when ${\displaystyle \psi =0}$. If ${\displaystyle \psi >0}$ and ${\displaystyle \theta >0}$, bistability will further require that that damage increases with load (${\displaystyle \dot{\overline{z}}(x)>0}$) which may be untrue early in the infection if ${\displaystyle \xi }$ increases very fast with load.

Relaxing the three assumptions altogether yields the following.

where ${\displaystyle {\delta }_y}$ is the parameter which determines the rate of defense consumption. We will assume in the following that ${\displaystyle F}$ increases with both ${\displaystyle x}$ and ${\displaystyle z}$. As in our first analysis, bistability will be possible if some of the terms in Equation S1-7 are negative.

Let us first consider the case where damage increases with load (${\displaystyle \dot{\overline{z}}>0}$). The only negative terms in Equation S1-7 are then ${\displaystyle \dot{\overline{z}}(x)\frac{\mathrm{\partial }G}{\mathrm{\partial }z}(\overline{y}(x),\overline{z}(x))}$ (as in Equation S1-6) which quantifies the negative impact of damage accumulation on defense production, and a term corresponding to the rate of defense consumption, ${\displaystyle -{\delta }_{y}x}$. Defense consumption is, therefore, clearly a mechanism which facilitates bistability. The activation of defense production upon damage detection will rather make it more difficult, because ${\displaystyle \dot{\overline{z}}(x)\frac{\mathrm{\partial }F}{\mathrm{\partial }z}(x,\overline{z}(x))}$ is positive.

In case damage decreases with load (${\displaystyle \dot{\overline{z}}<0}$) because repairing mechanisms are made more efficient when load increases, the negative terms in Equation S1-7 are ${\displaystyle \dot{\overline{z}}(x)\frac{\mathrm{\partial }F}{\mathrm{\partial }z}(x,\overline{z}(x))}$ and, again, ${\displaystyle -{\delta }_{y}x}$. This time, the activation of defense production upon damage detection will facilitate bistability.

In summary, as long as damage increases with pathogen load, bistability will be made possible by any mechanism that reduces defense when the infection progresses, be it defense consumption or lowered defense production caused by the accumulation of damage.

Appendix 1—figure 1

Download asset Open asset

A transient SPPL is a load which stabilizes for a long but finite period of time. This situation somehow resembles that of a chronic infection, but control is only transient and the infection will eventually enter a final phase where it reaches high loads and finally kills the host.

The progression of the infection slows down when a transient SPPL is approached because the trajectory evolves in the neighborhood of an equilibrium point; the fact that the infection eventually moves away from this point demonstrates that this equilibrium is unstable. Therefore, the SPPL cannot be transient unless an unstable equilibrium exists with ${\displaystyle \tilde{x}>0}$, which in turns requires that the system is bistable.

In such conditions, the pathogen load typically first peaks, then decreases until it reaches a plateau (a local minimum) and finally increases again while the disease proceeds towards its final phase. We computed the SPPL is such situations as the load at the local minimum. Computing the duration of chronicity is not as straightforward. We considered that control starts when load peaks, because it corresponds to the time at which immunity manages to curb bacterial proliferation. We further considered that control ends when load, after having decreased to the SPPL, reaches again the peak load. We estimated the duration of chronicity as the time that elapses between these two moments.

In most conditions, increasing ${\displaystyle x_0}$, decreasing ${\displaystyle y_0}$ or increasing ${\displaystyle z_0}$ increases the SPPL and shortens chronicity (Appendix 3—figure 1). In summary anything that strengthens immunity does lower the SPPL and, most of the time, lengthen chronicity. Exceptions to this rule exist, though: Appendix 3—figure 1 indeed illustrates cases (e.g. when ${\displaystyle \psi =0.5}$) where the duration of chronicity decreases when ${\displaystyle y_0}$ is increased. It should be noted, however, that this unexpected result may originate from the way we have defined the duration of control. We chose to start the period of transient chronicity at the peak that precedes the SPPL but other methods could be imagined, and it may be that some of the results we have obtained on this quantity depend on how it is precisely defined.

When the SPPL is transient, both its value and the duration of control vary with initial conditions, because they both depend on how close the trajectory will get from the unstable equilibrium. When stable, the SPPL corresponds to a stable equilibrium. It, therefore, lasts forever (which in more realistic terms means that the host will die from something else than the infection) and does not depend on initial conditions. This provides a method to experimentally distinguish the two situations: if the SPPL increases with the dose of pathogen a host is infected by, it must be transient. This could be used as an indirect way to demonstrate first that the system constituted by the interacting host and pathogens is bistable, and second that damage caused by the infection hinders the production of immune defense (see Bistability requires that damage hinders defense production).

Appendix 3—figure 1

Download asset Open asset

We propose that anything that strengthens the host immunity decreases the SPPL and, in most situations, lengthens the duration of control. We explored this rule in more details by running a elasticity analysis of the SPPL and the duration of chronicity with ${\displaystyle x_0={10}^{-6}}$ and ${\displaystyle y_0}$ and ${\displaystyle z_0}$ set at their homeostatic state (see Appendix 3—figure 2 and Appendix 3—figure 3). For each parameter, we numerically computed the derivative of the SPPL using a five point approximation. In the case of parameter ${\displaystyle u}$, we computed right derivatives, so that the condition ${\displaystyle u>v}$ holds in our computation. Similarly, we did not compute derivatives with respect to parameter ${\displaystyle v}$ alone, but rather varied both ${\displaystyle u}$ and ${\displaystyle v}$ at the same time, keeping ${\displaystyle u=v}$, so that ${\displaystyle F}$ always increases with ${\displaystyle x}$. Overall, we found that increasing the production of defense (by increasing γ or α) and making those defense more efficient (by increasing ${\displaystyle \delta }$) lower the SPPL. Increasing the rate at which damage is repaired (by increasing ${\displaystyle \xi }$) does the same (Appendix 3—figure 2).

What determines variation in the duration control is more complex. For example, when ${\displaystyle \psi =0.5}$ parameter combinations exists for which the duration of control increases when immunity causes more damage (i.e. when ${\displaystyle \eta }$ increases, see Appendix 3—figure 3). This is reminiscent to the situation we discussed in the previous section, and again we must warn that different ways of calculating duration of control are possible and that this paradoxical result might come from the method we have chosen. It should be noted however that the parameter range over which the duration of control is not negatively related to the value of SPPL is extremely limited.

Appendix 3—figure 2

Download asset Open asset

Appendix 3—figure 3

Download asset Open asset

Louie et al., 2016 proposed that tolerance could be assessed by measuring how fast the time to death decreases when the dose of pathogen is increased. We tested this idea by computing the slope of the relation of LT on dose and then running an elasticity analysis on the slope of the regression line, similar to that we ran on SPPL, PLUD, and HR. More precisely, we estimated slopes by varying the average dose ${\displaystyle x_0}$ from 10^-7 to 10^-5 (with 50 values sampled) and computing corresponding LT. We then adjusted a regression line to the values obtained. We then computed numerical derivatives of the slope of the regression line respective to each of 14 parameters of the model, using a method similar to that we utilized for other quantities.

As any other parameter measured on WHD or survival curves, the slope of the regression of LT on dose depends on both tolerance and resistance parameters (see Appendix 5—figure 1). For example, an increase in tolerance to damage ${\displaystyle z_d}$ brings the slope closer to zero. This is in line with the expectations of Louie et al., 2016: the more tolerant the host, the less LT should vary with dose. But an increase in ${\displaystyle \delta }$, the efficiency of defense, a parameter which should clearly be categorized as controlling resistance rather than tolerance, yields a steeper regression line.

Finally, increasing ${\displaystyle \xi }$, the rate at which damage is repaired, also makes the regression line steeper. Being able to repair damage efficiently should probably be considered as part of what makes a host tolerant to the infection. It is, therefore, surprising that ${\displaystyle \xi }$ and ${\displaystyle z_d}$ have opposite influences on the slope of the regression of LT on dose. We found that the influence of ${\displaystyle \xi }$ is lowered when ${\displaystyle \psi }$ is decreased (see Appendix 5—figure 1, ${\displaystyle \psi =0.5}$) which indicates that the unexpected effect of ${\displaystyle \xi }$ could stem from the negative impact of damage on the production of immune defense.

Overall, the slope of the regression of LT on dose cannot be considered as a simple proxy for tolerance. We found that it is sensitive to the complex interplay between damage accumulation and defense production, which makes it complicated to interpret.

Appendix 5—figure 1

Download asset Open asset

In any experiment of the type we described above, some deaths may occur for reasons other than infection. Measuring the PLUD on individuals that were not killed by the infection has of course no biological relevance. It may even cause spurious results if ever the risk of such death differs between treatments: treatments where individual tend to die before being killed by the infection should have lower PLUD, which may lead to flawed conclusions.

One way to limit this risk would be to remove extremely low PLUD values from experimental data. This can be done by fixing a somewhat arbitrary PLUD threshold; it could also be achieved more rigorously by using statistical tests specifically designed to detect outliers. In the following simulations, we evaluated the efficacy of this latter approach by using the Rosner test, which assumes that samples of size ${\displaystyle n}$ contain ${\displaystyle n-k}$ data points which are drawn from a Gaussian distribution while the ${\displaystyle k}$ remaining points are outliers.

For this purpose, we first used our model to simulate PLUD data. We then assumed that for each individual time to death is either determined by how fast the infection proceeds, as predicted by our model, or rather randomly drawn from a Weibull distribution, when death is caused by events that are independent from the infection. If the model predicts that death should occur later than this randomly drawn time value, we considered that the individual has not died from the infection. In such cases, we used our model to simulate the pathogen load at the randomly drawn death time.

Using a Weibull distribution allows us to consider cases where the background mortality rate does vary over time, during the course of experiment. We considered here a shape parameter of 1 (constant background mortality, yielding exponential distribution), 3 or 6 (meaning that mortality increases over time). Having fixed the shape parameter, we adjusted the scale parameter of the Weibull distribution so that we could control the proportion of individuals dying from causes other than infection. We varied this proportion from 2.5% to 25%.

We adjusted model parameters so that the distribution of PLUD values approximately matches the one we observed in our experiment (see Appendix 7—figure 2). We also considered that part of the variance in PLUD originates from uncertainties in load estimations. To reproduce this, we considered that actual load is 10⁷ and explicitely simulated the procedure of plating samples. More specifically, we used the dilution scheme presented above to predict the average number of CFU in a 5 μL droplet, and randomly drawn actual CFU numbers from a Poisson distribution. We then considered that only CFU numbers below 50 are countable (as in our experiments) and estimated load using a right-truncated Poisson model.

We found that the Rosner test adequately detects most of the individuals which died before being killed by the infection (Appendix 7—figure 2). Nevertheless, approximately 10% of them were not detected, whatever Weibull distribution we used. We evaluated how this lack of power could impact subsequent PLUD comparisons, by running simulations where a Wilcoxon rank test is used to compare two 50 individuals samples, one where all individuals died from the infection and another where some individuals died from other causes. All infections having the same dynamics in these simulations, the PLUD should not not differ among samples. But Appendix 7—figure 3 demonstrates that the presence of individuals who died before being killed by the infection creates spurious differences. As expected, increasing the proportion of such individuals worsen the problem. Conversely, increasing the shape parameter of the Weibull distribution tend to reduce the difference between samples because individuals tend to die late, even when not killed by the infection, at a load which, therefore, approaches the PLUD.

Appendix 7—figure 1

Download asset Open asset

Appendix 7—figure 1—source data 1

Protocol of dilution to estimate the Pathogen Load Upon Death (PLUD).

https://cdn.elifesciences.org/articles/104052/elife-104052-app7-fig1-data1-v2.xlsx

Download elife-104052-app7-fig1-data1-v2.xlsx

Appendix 7—figure 2

Download asset Open asset

Appendix 7—figure 3

Download asset Open asset

The second row of Appendix 7—figure 3 demonstrates that the problem is almost suppressed when outliers detected by a Rosner test are removed before the Wilcoxon test is run. Spurious significant differences can still occur with a 10% probability, but only when one sample contains more than 10% of the individuals which do not die from the infection. Using adequate controls in the experiment should allow to detect such situations: indeed, only individuals from the second sample should die in the absence of infection.