key: cord-0727821-cszogpdx authors: Shinbrot, Troy title: Immune dynamics in a time of covid date: 2021-06-02 journal: bioRxiv DOI: 10.1101/2021.06.01.446677 sha: 562cdd2c7444cbd89cdbcd7555b8c340f0892eb7 doc_id: 727821 cord_uid: cszogpdx Motivated by curiosities of disease progression seen in the coronavirus pandemic, we analyze a minimalist predator-prey model for the immune system (predator) competing against a pathogen (prey). We find that the mathematical model alone accounts for numerous paradoxical behaviors observed in this and other infections. These include why an exponentially growing pathogen requires an exposure threshold to take hold, how chronic and recurrent infections can arise, and what can allow very sick patients to recover, while healthier patients succumb. We also examine the distinct dynamical roles that specific, “innate,” and nonspecific, “adaptive,” immunity play, and we describe mathematical effects of infection history on prognosis. Finally, we briefly discuss predictions for some of the effects of timing and strengths of antibiotics or immunomodulatory agents. Coronavirus infections provide examples of several curious behaviors that highlight the remarkable dynamics of the immune system. Some of these can be explained through a basic understanding of how the immune system confronts an invading pathogen; others cannot. For example, it is not obvious why pathogens that grow exponentially only take hold if exposure exceeds a critical value 1 . Nor is it clear why seemingly similar patients with comparable viral loads produce widely diverging outcomes, e.g. one asymptomatic and another deathly ill 2,3 . Nor is it known why some patients steadily improve after initial infection, while others worsen, improve, and then relapse 4 , sometimes repeatedly 5, 6 . We examine a simple model that accounts for these and other enigmatic dynamical phenomena. We base our model on a straightforward predator-prey framework, in which immune agents (predators: ! ) attack pathogens (prey: ). Predator-prey models for immune response have been described before 7, 8, 9, 10 ; here we use these models to examine in immune dynamics in detail. We consider a version in which and ! change in time, , according to the schematic prescription: This in short says that pathogens reproduce, but are depleted by both adaptive and innate immune responses, and the immune system grows in response, but is self regulated by negative feedback. Most of these terms are straightforward to assign, as we describe next. Before we do so, we emphasize that what we present is only a simplified mathematical model: it defines the dynamics of pathogen growth competing against idealized innate and adaptive immunity, but entirely neglects real and important biological effects known to produce complicated disease behaviors. For example, interactions between immune function and inflammation 11 , as well as autoimmune disorders 12 , are topics in their own right. Many pathogens, moreover, have evolved devious mechanisms to battle the immune system. Tuberculosis, for example, defeats a key intracellular process so that it can hijack and reproduce within the very cells that destroy most invaders 13 . Helicobacter pylori constructs a robust biofilm to isolate itself from harsh environmental and immunological attack 14, 15 . The malaria parasite creates a vacuole on the membrane of liver cells, where it quietly reproduces insulated from the immune system so that it can periodically disgorge broods of new infective cells 16, 17 . And HIV relentlessly mutates over the course of years until it has amassed an army of different strains that overwhelm the immune system 18 . All of these effects produce complex disease dynamics, none of which are included in the present study. The value of the model that we present is not that it comprehensively describes the biology of pathogens and immune responses, but rather that it describes the mathematically necessary dynamics that must be present with or without other more elaborate biological effects. As we will see, these basic ingredients are sufficient to account for many curiosities seen in coronavirus and other infections. Pathogen kinetics: Pathogens typically reproduce at fixed intervals, so in Eq. [1] we set: where ! is a fixed growth rate. This produces growth defined by !" !" = ! , with exponential solution = ! · , where ! is the initial pathogen load. As for the remaining terms in Eq. [1], the immune system involves numerous interacting chemical and biological agents that are categorized in terms of whether they are adaptive (specific to a particular pathogen), or innate (nonspecific). For simplicity, we describe the adaptive response as a first order interaction, meaning a single immune cell destroys a single pathogen at a fixed rate, ! , and we describe the innate response to be zeroth order, meaning it depletes any invading pathogen at a fixed rate, ! : = ! . [5] This completes the definition of pathogen growth: As we have emphasized already, this is by no means a comprehensive description of realistic immune response, but as a simplified model we will see shortly that it captures the essential effects of these two immune mechanisms. For a more thorough description of the intricate workings of the immune system, an enjoyable overview can be found in Ref. [19] . As for pathogen growth described by Eq. [2], the adaptive immune system both trains cells to recognize new pathogens and recruits cells to attack known pathogens (the innate system does so as well, but to a more limited extent that we neglect here). Two central features characterize the adaptive response to pathogens: first, both training and recruitment take time, typically several days, and second both responses are "cooperative" in the biological sense, meaning that a mild infection provokes a mild response, but as the infective load grows, the immune response ramps up dramatically. Cooperative responses are often described by a "Hill" equation 20 , [7] where the exponent defines how cooperative the system is -i.e. how rapidly the response grows with increases in . A similar response has been used previously to model white blood cell production, which of course is central to adaptive immunity 20 . So we define the growth in immune response to obey: Here ! is a constant defining the rate of immune growth, and the ! term reflects the fact that immune cells themselves are required to produce that growth. We use = 2, which is most common in cooperative systems, but for the most part (an exception is mentioned at the end of our analysis) changes to affect our results little. The Hill term explicitly depends on the pathogen load delayed by the time needed to train or recruit a response. We would expect a recently encountered pathogen to produce a rapid response (small ); a novel pathogen would require a longer training time (large ). Finally, once a pathogen has been removed, the immune system returns to a quiescent state, so we include an autoregulatory term: [9] where the constant ! defines the speed of return to a number, ! , of adaptive immune cells. This completes the description of the immune response: [10] Armed with Equations [6] and [10] , we can now analyze how the model immune system can be expected to behave. We begin with the simplest cases and add complications sequentially. The simplest two cases are (1) a system lacking any immune function, followed by (2) innate immune response alone. We examine these cases by setting ! = ! = 0, which eliminates adaptive immune response, and for illustration we choose representative parameters ! = 2, ! = 1, ! = 0.5. This transforms Eq's [6] & [10] to: . [12] In this case, the pathogen and immune variables are independent of one another, and the equations can be solved separately, yielding two exponentials: where ! and !" are respectively initial values (numbers or concentrations) of pathogen and adaptive immune agents. It is useful to plot trajectories for these two first cases as shown in Fig. 1 Since !" = 0 along the nullcline defined by Eq. [15] , trajectories are entirely vertical along that curve (green in the figure), and similarly trajectories are entirely horizontal along nullcline [16] (red). All simulations in this work are performed in Mathematica, using its StreamPlot and NDSolve functions. From Fig. 1(a) , we can conclude, logically enough, that lacking all immune function any level of pathogen will grow -and from Eq. [13] , the growth will be exponential in time. Without loss of generality, we arbitrarily choose › 2 as a dangerous pathogen load, indicated by red shading in the figure. Adaptive immunity has no effect on the pathogen since we have set ! = ! = 0 in this example. Nevertheless trajectories approach a stable immune response at ! = ! , meaning that once interactions with pathogens have been eliminated, the last term in Eq. [10] does successfully autoregulate the concentration of immune cells so that it always reaches ! . There is no horizontal component on the green nullcline. Notice that in panel (a), with no immune response, pathogens grow for any initial or ! , and in panel (b) with only innate immunity, a critical pathogen load ! emerges (the green nullcline), below which the pathogen can be defeated. Initial pathogen loads within aqua region are cured; those within white region lead to unbounded pathogen growth. Intersections between the nullclines represent fixed points where neither nor ! change. Fig. 1(b) , on the other hand, shows two things. First, the intersection between the green and red nullclines indicated is a point where both = 0, so neither nor ! changes there. We call this a fixed point. Second, the green, = 0, nullcline separates pathogen loads that exponentially decrease (to the left) from those that exponentially increase (to the right). This means that the innate term produces a critical pathogen load, = / , below which the pathogen will die out, and above which the pathogen will grow. We will see this effect in more complicated examples in following sections. Already these simplest two examples provide an explanation for a first curiosity of immune function: even though infections typically grow exponentially, steady, nonspecific depletion of pathogens is sufficient to overcome minor infections. Moreover, the size of the infection that can be tolerated is determined by the ratio [15] of the rate of innate depletion, ! , to the rate of pathogen growth, ! . For this reason, we can tolerate a relatively large pathogen load of slowly growing pathogens, but rapidly growing pathogens reduce the critical pathogen load that can be destroyed. This is why we can tolerate limited exposure to Covid and other infections, whether we have previously been exposed to them or not. Eq. [15] also specifies that the limit to this exposure diminishes in a predictable way as the infection becomes more virulent. We turn next to examining how the adaptive part of Eq's [6] & [10] behaves, considering first the adaptive system by itself, and then considering the adaptive and innate systems combined. The innate system is defined by a single term, ! , in Eq. [6] . Setting that to zero and using representative values ! = 2, ! = 1 transforms Eq's [6] & [10] to: For the time being, we set the adaptive setpoint ! = 0.5; we also neglect the effect of the time delay, , needed to train and recruit adaptive cells: we will examine effects of these terms later. As before, it is useful to plot the nullclines, which provide a framework to understand the competing behaviors of pathogen and immune systems. Again, we do this by simply setting Eq's [17] and [18] equal to zero and solving: [20] In Fig. 2 , we plot these nullclines along with typical trajectories for two different values of ! , the autoregulatory rate at which adaptive immune cells return to the nominal setpoint, ! = ! . Notice that return to the setpoint can only fully occur once = 0, otherwise the first term in Eq. [18] will cause growth in ! . Notice also that for smaller initial adaptive immune response ( ! ‹ ! ), any infection will become worse before it becomes better. We will see that this is not always the case: some infections (e.g. those above the green nullcline, which have ample immune protection) only improve, while others improve and then worsen: we will introduce these later. From examination of Eq's [19] & [20] , one fixed point is = 0, ! = ! , indicated by an open black circle in Fig. 2 . A second fixed point is where the two nullclines intersect, as we saw in Fig. 1(b) ). We identify this fixed point by a solid circle, and behavior near this point plays an important role in determining outcomes of infections. In this case, notice in Fig. 2 (a) that the sample trajectory highlighted in blue spirals inward toward the fixed point. Because trajectories near the fixed point travel continually inward, the fixed point is stable: any disturbance will return the state to the fixed point. This fixed point occurs at a nonzero pathogen load, so the model apparently describes a chronic condition in which initial extremes of pathogen concentration diminish, but the pathogen is never fully destroyed. Thus the presence of chronic infections can also be understood using this model. At its heart, the chronic state arises because growth of the pathogen competes against depletion by adaptive immunity, but since that immunity diminishes with pathogen load, below a critical load the immune system becomes too weak to defeat the pathogen. We'll see shortly how the model nevertheless can describe eradication of a chronic infection. Additionally, because the blue trajectory in Fig. 2 (a) repeatedly spirals around the fixed point, it represents recurrence of increased and decreased pathogen load. Many diseases have recurrences with well-identified causes -for example tertian malaria is associated with a 48-hour rhythm produced by the reproductive cycle of the underlying parasitic infection 17 . Distinct from such rhythms caused by the pathogen, there are rhythms caused by the immune system: for example fevers commonly rise and fall daily in sync with circadian physiology. We emphasize that the recurrent spiral shown in Fig. 2 (a) is independent of rhythms imposed by pathogen, fever or normal physiological rhythms: it is intrinsic to the competition defined by Eq's [6] & [10] . So cyclic and recurrent infections are also expected outcomes from our model. In Fig's 2 we also display the effect of changes to immune autoregulation. Fig. 2(b) shows that an increase in the autoregulatory depletion rate, ! , of adaptive immune agents leaves qualitative behaviors of the pathogen-immune competition intact, but increases the chronic pathogen level -that is, the fixed point indicated by the solid circle moves to the right from Fig. 2 (a) to (b). Trajectories spiral around that fixed point as before, but along the way reach much higher pathogen concentrations. So if we again say that pathogen loads beyond the right side of the plots shown are life-threatening, then the blue trajectory at the slower immune depletion rate of Fig. 2 (a) would represent undesirable, but tolerable, disease levels, while the faster depletion in Fig. 2 (b) would, starting from the identical initial state, produce a dangerous infection. This means that rapidly suppressing the immune system, e.g. by abrupt pharmaceutical intervention, can lead to dramatic worsening in disease progression. The root cause of this is that our model dynamics have an intrinsic time scale, and holding the pathogen in check relies on the enduring presence of immune agents. Reducing their effect too rapidly can allow the pathogen to grow too fast to be brought under control. Contrariwise, Fig. 2 (c) shows that lowering the rate, ! , of immune autoregulation can be expected to reduce illness and lower the chronic pathogen load. In that figure, we plot trajectories for the extreme case without any autoregulation. In that case, all trajectories end where they hit the vertical axis, i.e. at = 0, ! › ! . Immune agents remain at this unregulated level forever, and so the organism must tolerate an elevated immune level for the rest of its existence (a metabolic expense that evolution has presumably optimized). To recap, immune autoregulation in our model prescribes that adaptive immune levels will tend to return at rate ! to the setpoint ! . In Fig. 2 we considered how the rate ! affects the system; in Fig. 3 (a) & (b), we overview effects of the setpoint. In Fig. 3 (a) we plot trajectories with the immune setpoint reduced, to ! = 0, -i.e. where immune cells can respond to an infection, but are afterward autoregulated to a state with zero residual cells. This case is instructive because it is conservative, which implies that trajectories always travel in closed orbits. As an example, the magenta trajectory in Fig. 3 (a) will periodically cycle around the fixed point = 1, ! = 2 forever. In fact almost all trajectories, for any initial state, encircle that fixed point. The blue trajectory, starting at low pathogen load and finite immune response, also does so, but it produces very large pathogen loads on its way in an extensive voyage around the fixed point. Thus just as in Fig. 2(b) , where immune cells were rapidly depleted, reducing the setpoint suppresses the immune response and allows pathogen levels to grow without significant opposition. Fig. 3(b) , however, shows that unlike an increase in the rate of immune depletion, an increase in the setpoint strengthens the stabilization of the fixed point, causing the immune response to grow faster than the pathogen. As a consequence, starting from the same initial state as in Fig. 3(a) , rather than producing a dangerous infective growth, the blue trajectory now shows a rapid approach to the chronic fixed point. , we see that increasing the immune setpoint provides substantially better protection against infection than changing the autoregulatory rate. Apparently, and logically, the presence of significant residual levels of immune agents is highly beneficial to the function of the adaptive immune system. This is the mathematical interpretation of the value of prior exposure or of vaccination. Increasing the rate, ! , of autoregulatory control over homeostatic immune levels can reduce the metabolic expense of maintaining these agents, but without this expense, the agents will senesce or otherwise be cleared, allowing an incipient infection to grow to dangerous levels before the immune system has time to rally a response. We have looked so far in this section at autoregulatory control over immune response. This is encapsulated in the last term in Eq. [10] , which depends on ! and ! , but plainly the adaptive immune growth, specified by ! , competes against this autoregulatory depletion. To explore the effect of this competition, in Fig. 3 (c) we increase that rate from its nominal value ! = 1 up to a new value, ! = 4: twice the rate of pathogen growth, ! = 2. The other parameters are returned to their nominal values originally shown in Fig. 2 (a): ! = 2, ! = 0.5, ! = 0.5. As Fig. 3(c) shows, by increasing ! we achieve much the same end as increasing the immune setpoint, ! . This makes sense: one can achieve a similar goal either by starting at a higher immune level or by driving the system to that level faster. Fig. 3 (c) also shows, however, that faster immune growth causes substantial overshoot, producing large excursions in immune levels, and consequently the system takes longer to reach steady state. Summarizing findings up to this point, we have seen that the innate immune system by itself can overcome minor infections. Mathematically speaking, as it does so an unstable demarcation (the green nullcline in Fig. 1(b) ) emerges between tolerable and dangerous infections. The adaptive system alone can constrain larger infections, but it tends to stabilize those infections, leading to a chronic state (that may or may not oscillate depending on autoregulatory details). The interplay between unstable and stable states is central to understanding the combination of innate and adaptive immunity, and in Fig. 4 we show several different behaviors produced by this interplay as the innate strength, ! , is increased. Under identical conditions but increasing the innate immunity coefficient to ! = 0.2 makes the fixed point unstable. As a consequence, the same initial condition that was attracted to the fixed point in panel (a) is now repelled, leading to zero pathogen load as shown in solid blue trajectory. As in prior plots, infections that start within the aqua region lead to cure; infections within the yellow region do so as well, but the disease gets worse before it gets better. Very low immunity, in white, gets steadily worse. (c) Modeling a boost to the immune response by further increasing ! to 0.5 causes the green nullcline to move to higher pathogen loads. This nullcline divides regions that steadily improve from those that worsen before they improve. So increasing innate immune function -i.e. increasing ! -causes the same initial state as shown in solid blue in panels (a) and (b) to be rapidly cured. Starting with Fig. 4(a) , we show trajectories produced with ! = 0, i.e. without innate immune function. As in the earlier Fig. 2(a) , we use illustrative values ! = 2, ! = ! = 1, = 0, ! = ! = 0.5 in Eq's [6] and [10] . In Fig. 4(a) , we identify the stable fixed point that attracts nearby trajectories as a solid magenta circle, and we show a representative orbit around the fixed point in blue. Beginning at the same initial point, Fig. 4(b) shows the solid blue trajectory in the presence of both adaptive and innate immunity, produced by setting ! = 0.2. Notice that this trajectory still encircles the fixed point, but rather than being attracted, it is repelled by that point -hence we now identify it with an open circle using the convention introduced in Fig. 2 . That is, introducing innate immunity has caused the fixed point to transition from stable (attractive) to unstable (repulsive). This occurs in this example when the innate immunity strength in Eq. [6] , ! , exceeds about 0.1018. When the fixed point becomes unstable, it pushes the solid blue trajectory away as we see in Fig. 4(b) , so that it now reaches = 0 -which is to say so that all pathogens are destroyed. This means that combining innate and adaptive immunity allows more serious infections to be controlled, though neither immune response alone can do the job. The innate system alone can eliminate minor infections but cannot control larger ones. The adaptive system alone can stabilize larger infections, but tends to reduce its intensity as the infection diminishes -and hence can produce chronic conditions. We will see that the delay of the adaptive system described by Eq. [10] is also important to chronic conditions, but for now we settle for the observation that combining innate and adaptive responses can effectively combat infections that neither can cure alone. Notice that in Fig. 4(b) there are two qualitative outcomes. First, trajectories beginning in the shaded aqua region above and left of the green nullcline lead monotonically to decreased pathogen load -that is, outcomes steadily improve. An example is shown as a dashed blue line. Trajectories in this aqua region that begin to the left of the red nullcline see a decrease in both immune and pathogen levels; those above the green but to the right of the red nullcline require an increase in immune expression to destroy the pathogen. In either case the pathogen is steadily and completely removed. Second, trajectories like the solid blue line that start in the yellow region produce an increased pathogen level before the pathogen is removed, meaning the disease gets worse before it gets better. We have seen this before in Fig's 2 and 3 -we will see shortly, though, that some diseases in this model can do the opposite: i.e. get better and then get worse. Before we examine more complicated cases, note that even stronger innate immunity drives green nullcline further to the right: this is shown in Fig. 4(b) , where we raise ! from 0.2 to 0.5. To stretch the point, this could also be considered to result from an adjuvant that boosts a weak immune response. The green nullcline, recall, distinguishes trajectories that monotonically improve from those that first worsen: this occurred as well in Fig. 1(b) without adaptive immunity. A companion finding with combined innate and adaptive community is that the green nullcline can move past an initial state and so switch the outcome between one that steadily improves and one that worsens. For example the solid blue trajectory in all panels of Fig. 4 starts from the same location, but leads in panel (a) to a recurrent infection, in panel (b) to an infection that worsens and then resolves, and in panel (c) to an infection that rapidly and monotonically improves. Thus a moderate increase to model innate immunity can resolve chronic infection following a period of recurrence: a larger increase can rapidly resolve the same infection. We showed in Fig. 4 that adjacent initial states can lead to steadily worsening (white), improving (aqua), as well as chronic or recurrent (sickly green) infections. We also mentioned that some states worsen and then improve, while others improve and then worsen. In Fig. 5(a) , we identify regions, or "basins" that lead to each of these outcomes for the combined innate and adaptive immunity case shown in Fig. 4(b) . Fig. 5(b) shows an enlargement near the central fixed point: states within the black basin ultimately reach the redline indicating a dangerous pathogen load (i.e. patient death). In both panels, colored regions ultimately get better ( → 0); colorless regions ultimately get worse › 2 . The boundaries between distinct basins interleave closely, so small differences in initial state can produce large differences in outcome. Fig. 4(b) , color coded based on outcome. At the highest ! values (aqua), essentially any pathogen load produces steady improvement. In the lower ! band identified in black, the disease improves ( decreases) and then worsens ( increases). Lower ! still (yellow) produces worsening then improvement, and the lowest ! (gray or white) produces steady worsening. (b) Enlargement of region near fixed point shows that the outcome can ultimately worsen or improve depending on exact initial state. Initial conditions (black) that ultimately worsen spiral inward toward the fixed point. Initial conditions within colored regions ultimately improve. The strong dependence of outcome on initial state can be traced to the root cause, shown first in Fig. 2 , that there are two fixed points in our simplified model. We highlighted the importance of stability of a first fixed point in Fig's 4 and 5 ; the second fixed point is circled in Fig. 6 closer to the origin. Two nearby trajectories that approach this fixed point are shown in blue in Fig's 6(a) and (b) ; both panels use identical parameters as used in Fig. 4(b) . A first trajectory, in Fig. 6(a) , leads to a rapid elimination of all pathogens, while a second, in Fig. 6(b) , produces an exponential growth in pathogen levels. The parameters are exactly the same in both panels, but pathogen loads initially differ here by a part in 1000. This deviation can be arbitrarily small, for instance initial values of , ! = 0.3824534 , 2 and 0.3824535 , 2 differ by about 3 parts in 10 million (approaching the limits to computational precision), and again the first leads to a rapid cure while the second leads to dangerous pathogen growth. It is not difficult to show that the circled fixed point is always unstable 22 in this system, and this sensitivity to initial infection load is intrinsic to systems with unstable fixed points such as this. Thus the finding that very similar patients can experience very different outcomes is a natural expectation of the simplified model. , ! = 0.3826,2 , which shows exponential pathogen growth. In both panels, the same parameters are used as in Fig. 4 We learned from Fig. 4(b) -(c) that changing parameters such as rates of growth or regulation can move nullclines and so alter a disease's progression. From Fig. 5 , we saw that the same can be accomplished by changing the state, i.e. the pathogen load, , or immune response, ! . This raises the obvious question of whether outcomes can be intentionally controlled by making changes to parameters or variables. For example, effects of antibiotics can be modeled by recalling that the final, innate, term in Eq. [6] reduces infections at a fixed rate, and mathematically speaking this could just as well be caused by exogenous drug treatment as by endogenous immunity. As a case in point, it is clinically recognized that early intervention is important in the treatment of many diseases: a fact that is commonly ascribed to causes including the reduction in inflammation, tissue damage and biofilm formation by the offending organism. In addition to these very real causes, a fundamental dynamical effect is at work, as illustrated in Fig. 7 . In both panels, we plot the movement of the !" !" nullcline produced by an increase in ! from 0.1 to 1 -representing a faster depletion of a pathogen that might be produced by an antibiotic treatment. In panel (a) we illustrate an early intervention, showing treatment when the state has reached the red star, to the left of the moved nullcline. In panel (b) we show the same treatment applied later, when the state has reached the red star at higher pathogen load. The resulting trajectories in both cases are shown in blue, and follow the underlying trajectories, both at ! = 0.1 and ! = 1 shown in gray. Evidently, early intervention can produce a cure that later intervention -even using a stronger treatment -may fail to achieve. An essential feature of adaptive immunity is that it takes time to train or recruit a specific immune response 20 . The delay, , is included in Eq. [10] to account for this time, however including carries a mathematical burden: in order to solve for the state at a time, , one must know the history of ( ) and ! ( ) at all times between − and . Since there is an infinite number of such times, such "delay differential equations" are formally infinite dimensional, and consequently it is not possible to provide a uniformly applicable solution. In practice, the extent to which delay differential equations depend on history varies considerably, and the study of these equations is a topic unto itself 23 . To make numerical solutions tangible, we will consider two of the infinite number of possibilities. First, we will examine the situation in which a disease initially changes slowly compared with the time delay -that is, we define the initial pathogen history to be constant, = ! for -‹ ‹ 0. This is a simple case that can be reasonably cleanly analyzed. Second, we will examine a pathogen that grows exponentially at the fixed rate ! up to the value ! , i.e. = ! ! ! ! during the time -‹ ‹ 0. This case is more realistic, but since the pathogen load changes in time before the immune system has been engaged, it produces hysteretic effects that can seem paradoxical. We start by considering the case shown in Fig. 4(c) , with strong innate immunity ( ! = 0.5), and the other parameters as before: ! = 2, ! = ! = 1, = 0, ! = ! = 0.5. Without no delay, = 0, we reproduce the earlier results in Fig. 8(a) and superimpose a new representative trajectory in blue. The initial state shown produces a worsening of infection followed by a reduction in both and ! . This is a satisfactory outcome: the pathogen is eliminated, and autoregulation reduces the immune response to close to its original level. As the delay is increased, the immune response is further retarded, which allows the pathogen to grow for a time before the adaptive system responds. As a result, the maximum pathogen load, !"# in Fig. 8 , increases with delay time. Fig's 8(b) & (c) show trajectories starting from the same initial state as in Fig. 8(a) , but with delays = 1 and = 2 respectively. At = 1 (Fig. 8(b) ), !"# increases modestly; more significantly the adaptive immune response begins later in the game, and so the ultimate value of ! at which the pathogen is destroyed is much larger than before. This is beneficial from the standpoint that the immune system ends up highly energized and ready to combat a similar infection, but also potentially results in over-activation of the immune system. To make clear what the immune-pathogen competition is doing both in the "pre-history" stage, before the simulation begins at = 0 and afterward, we color the trajectory cyan during the period when the initial pathogen is held constant (-‹ ‹ 0), and blue thereafter › 0 . The cyan part of the trajectory depends on the prior history; the subsequent trajectory in blue is reliable given the history shown in cyan. We will describe some aspects of how history can affect outcome shortly. At longer delay yet, = 2 in Fig. 8(c) , !"# grows beyond the redline, indicating a dangerous pathogen load. This is to be expected: essentially the immune system is retarded beyond the time when it is needed, so can grow almost unconstrained. This fact allows us to redefine the basin of attraction for states that can be cured as a function of delay. So for example, recall from earlier figures (e.g. Fig. 4 ) that we identified initial states in yellow and aqua that lead to a cure ( = 0). As a delay is introduced, the curable region contracts as shown in Fig. 9 for delays ranging from = 0 to = 4. In Fig. 9 (a), we show the example of strong innate immunity, ! = 0.5, earlier considered in Fig. 4(c) . Initial states in the cross-hatched area would be curable with no delay, but as delay is introduced, the curable area becomes eroded -meaning that fewer states would be curable. Fig. 4(c) ), longer delays erode the cure basin: the cross-hatched region shows initial points that would be cured for = 0, but not for longer delays. Notice that a cusp grows as the delay increases, and so for nonzero delays, an immune response such as shown by the dashed line results, moving left to right, in a cure for small , death for larger , cure for still larger , and death for the largest values. (b) Weaker immunity: Using ! = 0.2 (cf. Fig. 4(b) ), the region of initial points that are cured shrinks. Note that with no (or small) delay, nearby initial states that lead to cure or death interweave as shown in Fig. 4 , but is eliminated by either stronger immunity (panel (a)) or significant delay (this panel). (c) Sample trajectories: Four trajectories starting from initial points indicated, using ! = 0.2, = 2, confirm that states evolve sequentially to 1: cure, 2: death; 3: cure, and 4: death. Notice as described in text that the delay term in Eq. [10] implies that trajectories with different histories can cross one another, which does not occur without delay, as in Fig. 4. As one might expect, when adaptive immunity is weakened, the range of states that can be cured is diminished, and even more so with delayed response. This is shown in Fig. 9(b) , where we plot basin boundaries as in 9(a), but with immunity reduced from a stronger value, ! = 0.5 to a weaker one, ! = 0.2. The cross-hatched area grows, indicating that a larger range of states is sensitive to delay when the immune system is weakened. More interestingly, a cusp develops as identified in Fig. 9 (a). Cusps and their consequent "catastrophes" have been well studied 24 ; the result here is that for a fixed immune response -for example, ! = 1.35 as indicated by the dashed line in Fig. 9(b) , increasing the pathogen load causes the ultimate outcome to switch from cure (at low ) to death (above about = 0.3) to cure again (~0.8) and finally to death (above ~1.6). The same occurs at stronger immunity ( Fig. 9(a) ), but the effect is less pronounced. For any immunity strength, this is a highly counterintuitive finding: it means that although very small or very large pathogen loads logically enough lead to patient recovery or death respectively, between those extremes one patient with a lower load may die, while an identical patient with a higher load may survive. The cause of this paradoxical result is that a higher pathogen load can provoke a stronger immune response that, combined with delay, can lead to a successful outcome, while a lower load can permit the pathogen to gain in strength, "sneaking up" on the immune system so to speak, and ultimately overcoming it. This may be made more evident in Fig. 9(c) , where we plot trajectories starting from four initial pathogen loads, : = 0.2, : = 0.7, : = 1.2, and : = 1.7, all for the same immune response, ! = 1.35. The first state is promptly cured, the second produces steady pathogen growth, the third worsens and then steadily improves, and the last rapidly dies. Notice that the trajectory labeled 2, in blue, crosses trajectories 3 and 4that is, starting from state 3 (with = 1.2 and ! = 1.35 held constant for time = 2) leads to a cure. This occurs because Eq. [10] : [21] holds !" ! constant, at a relatively high value, for the period from time = − until time = 0 when the simulation begins. Meaning that ! can grow at this constant value as soon as the simulation begins: as the trajectory beginning at state 3 shows. On the other hand, beginning at state 2 and passing through state 3, the simulation is held at the relatively lower value of that was established at state 2. Indeed, state 2 starts very near to the red !" ! !" nullcline, so ! changes by very little initially: as the trajectory beginning at state 2 shows. Thus the pathogen can grow, or sneak up on, the immune system, leaving it uninfluenced by the pathogen growth. We remark that it may appear from Fig. 9 (b) that the switching of cure and death outcomes is superficially similar either with or without delay. That is, if one travels along the dashed line in Fig 9( b) either = 0 or with = 2, one will in both cases encounter alternating regions leading to cure and to death. Mathematically, however, despite the apparent qualitative similarity, these situations are different. Without delay, trajectories are strictly two-dimensional, because the dynamical system defined by the ordinary differential equations [6] and [10] is itself two-dimensional, and resulting trajectories cannot cross one another. As we have remarked, however, the dynamics in the presence of delay is higher dimensional, and so trajectories in basins that lead to cure or death can cross one another (or more precisely, appear to cross in the two-dimensional projection of the underlying higher dimensional system). Crossings shown in panel (c) carry the meaning that given a state, , ! , at a time , different outcomes can emerge depending on the history prior to . Apropos of this observation, the specific outcome shown in Fig 9 obtains in our first, simpler, choice of history in which is held constant for time before beginning the simulation. A different choice of history would produce a different specific outcome. Nevertheless, it is a general result independent of choice of history that delays allow trajectories to cross so that a smaller initial pathogen load can produce a larger ultimate infection. To see this, we turn to the more realistic case of a pathogen that grows exponentially in time with rate ! , so that = ! ! ! ! during the time -‹ ‹ 0. Not surprisingly, this more elaborate history results in more elaborate trajectories. In Fig. 10 , we plot results using this exponential history with the same parameters as in Fig. 9(b) . We see that the relatively modest cusps produced previously become embellished into a spiral wave shape. In detail, as delay is increased, the original spiral shown in Fig. 4 without delay (gray in Fig. 10(a) ) unwinds, producing a basin boundary that bends through radians as shown at around = 1 (blue), about /2 radians near = 2 (red), and becomes single valued around = 4 (orange). That is, for large delays, corresponding to each adaptive response value there is a single critical pathogen load, ! , above which the pathogen will grow. The emergence of a critical load is similar to the case without adaptive response, except that the value of ! grows with ! as shown in Fig. 10(a) , and indeed the longer the delay, the more vertical (i.e. independent of ! ) the critical load becomes. This makes sense: as the saying goes, "justice delayed is justice denied," and in the same way delaying an adaptive immune response long enough is equivalent to removing it altogether. In Fig. 10(b) , we highlight the spiral wave shape for a delay = 1: evidently in this case, at a constant initial immune response ! = 2.2, the exponential growth of pathogen produces a sequence of five alternating outcomes. For in the three ranges in aqua ([0, 0.43] or [0.68, 0.89] or [0.98, 2]) a cure results; for between these ranges, in white, death results. Fig. 10 (c) displays trajectories, which again appear to cross, in each of these ranges. Evidently, both simplistic (Fig. 9 ) and more realistic (Fig. 10) history models produce paradoxical outcomes, so that for the same immune state a worse infection can produce a better result. We mentioned earlier that delay differential equations can be very complicated; a few examples of this are shown in Fig. 11 . First, we know from Fig. 4 that increasing the innate immunity (controlled by ! ) causes the stable fixed point to become unstable, and without a delay, we saw that this allows combined innate and adaptive immune systems to cure chronic diseases. It turns out that in the presence of delay 25 , trajectories become bounded rather than exiting the domain either below ! = 0, which we designated a cure, or above ! = 2, which we associated with death. We demonstrate that trajectories become bounded in the dark blue trajectory in Fig. 11 (a), which migrates inward as shown. In that panel, we use ! = 2, ! = 1.1, ! = 0.015, ! = 1, ! = ! = 0.5 and = 0.155. Using a different initial state, we get an expanding trajectory, shown as light blue in Fig. 11(a) . Since trajectories outside of the magenta orbit migrate inward, and trajectories inside migrate outward, it is a wellknown consequence of the Poincaré-Bendixson theorem that the magenta orbit itself must be a limit cycle. That is, for long times, all trajectories within a bounded region outside of the dark blue trajectory must approach the magenta orbit -and correspondingly, all states must periodically cycle indefinitely. Technically, the guarantee of a limit cycle assumes that the system is two-dimensional; computationally this is close enough to being the case that a limit cycle does in fact appear. For the same parameters as in Fig. 11(a) , but using smaller delays than about = 0.1, the central fixed point becomes stable, and a constant (rather than periodic) chronic state is approached. For larger delays the fixed point becomes globally unstable, and either a cure or death is reached as in previous examples. We note in both Fig. 11(a) and 11(b) , shown next, actual trajectories deviate from the streamlines shown in the background. This is a signature of the delay, which causes trajectories to follow derivatives obtained earlier in time than those produced by instantaneous evaluation. Other states are also easily produced; for example in Fig. 11(b) we show the opposite of a limit cycle: an unstable cycle obtained at ! = 3, ! = 1.1, ! = 0.08, ! = 1, ! = 0.5, ! = 1.325 and = 0.1. Trajectories beginning outside of that cycle expand to either reach a cure or death, and trajectories inside of the cycle orbit repeatedly and slowly approach a stable chronic state. It is unclear if this type of behavior is seen physiologically: it would represent a chronic state that, if significantly perturbed by treatment or environment, could either degenerate into a cure or death: a worrying possibility! Finally, we remark that other behaviors than those examined here are also possible. A couple of relatively simple examples that do not rely on delay can be generated by increasing the cooperativity, , from 2 to 4. This makes the red nullcline strongly sigmoidal as shown in Fig's 11(c) & (d) . The nullcline is weakly sigmoidal with = 2, but = 4 is preferable for illustrative purposes. We can lower the sigmoid on the ! axis by changing ! to 2 and ! to 1.2 but keeping the other parameters as before ( ! = 3, ! = 1.1, ! = 0.08, ! = 1 and = 0.1). This causes the two nullclines to separate as shown in Fig. 11(c) . This separation, referred to as a tangent bifurcation, eliminates intersections between the nullclines, and so the fixed points at the intersections vanish, and so does the possibility of chronic or recurrent states. Rather, as the blue trajectories starting from four different initial states illustrate, all trajectories end in death. Biologically, this occurs because the lower immune setpoint ( ! → 1.2 instead of 1.325) combined with the faster regulation to that setpoint ( ! → 2 instead of 0.5) causes the adaptive immune response to rapidly diminish as pathogen load increases. As shown in Fig. 11(c) , this forces all trajectories to quickly approach the saturation value of the ! nullcline ( ! = 1.38), irrespective of initial state. The tangent bifurcation shown in Fig. 11 (c) arose when the sigmoidal red curve was moved downward (by changing ! and ! at increased ). We can instead move that nullcline upward to produce a second tangent bifurcation, as shown in Fig. 11(d) . To do this, we keep = 4 and set ! = 1, ! = 3, ! = 2.75. The other parameters remain as before ( ! = 1.1, ! = 0.08 and = 0.1). As in Fig. 11 (c), these parameters are chosen for illustrative clarity, and are not special: many other values will also produce a sigmoidal ! nullcline separated above from the green nullcline. Whenever this separation occurs, all initial states end in a cure, as shown by the blue trajectories in the figure. This isn't mysterious: biologically, it arises because the pathogen growth rate, ! , was decreased from 3 to 1, and immune values are consequently forced rapidly ( ! → 3 instead of 2) to a high setpoint ( ! → 2.75 instead of 1.2). This represents a slowly growing pathogen attacked by immune agents that rapidly grow to a high value. Naturally enough, such infections are invariably cured. As we mentioned, the bifurcations shown in Fig. 11 (c) & (d) do not rely on delay; nevertheless they hint at other possible behaviors. The Mackey-Glass 20 equation defines a single variable that describes growth of blood cells, and when delay is included this produces chaotic and other complicated solutions. We have analyzed straightforward dynamics of two variables in our equations [6] and [10] , and these dynamics already exhibit plentiful dynamical complexities, many of which duplicate behaviors known to be present in immune responses, and a few of which may represent new unreported states. It seems inevitable that more detailed analysis including delay is bound, like the simpler Mackey-Glass system on which our model is based, to exhibit even more complicated solutions. To Physicists, it is surprising that infections that grow exponentially require an exposure threshold to persist, and that once an infection that could be successfully treated by antibiotics exceeds a larger threshold, no amount of antibiotic will prevent its growth. To Biologists, it is well known, but difficult to understand, why some infections persist, either in a relatively stable chronic form, or in a periodic state that regularly grows and shrinks. Likewise it is unexplained why some patients get worse before getting better, while others get better before getting worse. Nor why one patient can encounter a disease and survive, while a medically indistinguishable patient may succumb. All of these observations are explained by a minimal mathematical model of competition between a pathogen and an immune system consisting of both innate and adaptive parts. Specifically, we have found the following. First, considering the case where delay in mobilizing an adaptive immune response is neglected: 1) The innate immune system alone produces a critical pathogen load, ! below which the pathogen will die out, and above which the pathogen will grow. This provides an explanation for why we can tolerate limited exposure to infections. 2) The adaptive system alone produces an attractive fixed point with constant pathogen load and immune response. This can be interpreted as a chronic state. Moreover, autoregulation of the adaptive system can lead to recurrent cycles of infection. 3) Slow immunosuppression can reinforce these recurrent cycles, while fast immunosuppression enlarges the growth and shrinkage of pathogen levels during these cycles, leading potentially to dramatic pathogen growth. 4) Combining innate and adaptive immune responses allows more serious infections (that neither system alone could control) to be resolved. In short, the innate system can clear a minor infection, but cannot cope with a larger infection. The adaptive system can stabilize a larger infection, but cannot by itself destroy the infection. Combined, the two systems can both stabilize and eliminate larger infections. 5) Infections can exhibit several distinct types of dynamics. They can monotonically improve or worsen. They can worsen and then improve. They can improve and then worsen. And they can cycle repeatedly before ultimately resulting in a cure or in death. 6) The states, defined by initial pathogen load and adaptive immune status, leading to ultimate cure or death can be very closely interleaved, especially close to the chronic steady state. This has two implications. First, small differences in initial state can produce large differences in outcome. Second, in principle intentional intervention could steer the competition between pathogen and immune system to a state leading to a successful outcome. 7) With respect to interventions: a. Early intervention (e.g. with antibiotics) can produce a cure that later intervention may not. b. Indeed, once a pathogen level has exceeded a point of no return, no amount of intervention may effect a cure. c. As one would expect, vaccination or prior infection greatly improves the response to a new pathogen exposure. Second, if delay of the adaptive system is included: 8) Paradoxical outcomes arise, for example a patient subjected to a lower pathogen load may succumb while a comparable patient subjected to a larger load may survive. This appears to be due to the ability of a weaker pathogen load to grow to a dangerous level in the presence of a weaker immune system. This permits the pathogen to sneak up on the immune system, before the immune system has mobilized an effective response. The nonlinearity of the adaptive system produces a significantly stronger response to higher pathogen loads, which can produce the paradoxical outcome in which a sicker patient may survive while a less sick one expires. 9) These paradoxical responses appear to be generic in that simpler as well as more realistic delay models can both generate alternating cure and death outcomes as the pathogen load is increased. 10) Other, more complicated states are also possible, including stable limit cycles and chronic states that can degenerate into a cure or death. These are purely mathematical findings that follow from a straightforward model that utilizes a minimum of assumptions about the detailed workings of the immune system. We stress that our mathematical results can absolutely not be taken to suggest that more detailed and realistic workings of the immune system and of pathogenic evasion strategies are not important. Rather, our findings show that mathematics, on its own and without those details, provides mechanisms for a number of veridical dynamical behaviors that may otherwise be difficult to rationalize. Other possible complicated dynamics are also present in this system, including stable limit cycles and chronic states that can degenerate into a cure or death. Based on prior work of related equations, more complicated behaviors also seem likely. This work is partially supported by NSF CBET, award no. 1804286. The data that supports the findings of this study are available within the article. Experimental evidence of a pathogen invasion threshold Recurrence or Relapse of COVID-19 in Older Patients: A Description of Three Cases Repeated COVID-19 relapse during post-discharge surveillance with viral shedding lasting for 67 days in a recovered patient infected with SARS-CoV-2 Dynamics of recurrent viral infection Microbial dose response modeling: past, present, and future Predator-prey equations simulating an immune response Mathematical model of immune processes A predator-prey-disease model with immune response in infected prey The trinity of COVID-19: immunity, inflammation and intervention The autoimmune diseases How the immune system works, 4 th Edition Biofilm formation by Helicobacter pylori Biofilm, pathogenesis and prevention-a journey to break the wall: a review The parasitophorous vacuole of the blood-stage malaria parasite The malaria parasite has an intrinsic clock How HIV defeats the immune system How the immune system works, 4 th Edition Pathological conditions resulting from instabilities in physiological control systems Nullclines and phaseplanes Nonlinear dynamics and chaos with student solutions manual: With applications to physics, biology, chemistry, and engineering An introduction to delay differential equations with applications to the life sciences Catastrophe theory Predator-prey models with delay and prey harvesting Yong, E.(2020) Immunity is where intuition goes to die, The Atlantic https://www.theatlantic.com/health/archive/2020/08/covid-19-immunity-is-the-pandemics-centralmystery/614956/