\begin{matrix}c \\ ◯ \end{matrix} Revista de Matemática: Teoría y Aplicaciones 2020 27(1) : 141–156 cimpa – ucr issn: 1409-2433 (Print), 2215-3373 (Online)

doi: https:// doi. org/ 10.15517/ rmta. v27i1. 39965

backward bifurcation in neutrophil-pathogen interaction

bifurcación hacia atrás en interacción neutrófilo-patógeno

Christopher M. Kribs^{∗}  Omomayowa Olawoyin^{†}

Received: 12/May/2019; Revised: 31/Aug/2019; Accepted: 17/Sep/2019

∗University of Texas at Arlington, Department of Mathematics, Arlington TX 76012, United States. E-Mail: kribs@uta.edu

†Misma dirección que/Same address as: C. Kribs. E-Mail: omomayowa.olawoyin@mavs.uta.edu

142 c.m. kribs — o. olawoyin

Abstract

Bacterial infections elicit immune responses including neutrophils, whose recruitment is stimulated by the bacteria’s presence but which die after eliminating those bacteria. This dual interaction between bacteria and neutrophil concentrations, more complicated than the simple predator-prey relationship that describes macrophage-bacteria interactions, creates an environment in which neutrophils may only be able to clear sufficiently small infections. This study describes this relationship using a simple nonlinear dynamical system which exhibits bistability behavior known as a backward bifurcation. Bacterial growth is assumed limited by a key nutrient. In contrast to a previous study which held neutrophil and nutrient levels constant and required saturation terms to produce bistability, our model shows that simple bilinear terms support bistability when nutrient and neutrophil densities are allowed to vary in response to bacterial density. An example application involving Borrelia burgdorferi, which feeds on manganese, illustrates why neutrophils’ rapid response is key to their ability to contain bacterial infections.

Keywords: bistability; resource limitations; phagocytosis; hysteresis.

Resumen

Las infecciones bacterianas provocan respuestas inmunitarias, incluyendo neutrófilos, cuyo reclutamiento es estimulado por la presencia de la bacteria pero que muere después de eliminar esas bacterias. Esta doble interacción entre las concentraciones de bacterias y neutrófilos, más complicada que la simple relación depredador-presa que describe las interacciones entre bacterias y macrófagos, crea un ambiente en el que los neutrófilos tal vez sólo puedan despejar infecciones suficientemente pequeñas. Este estudio describe esta relación utilizando un sistema dinámico no lineal sencillo que exhibe un comportamiento de biestabilidad conocido como una bifurcación hacia atrás. El crecimiento bacteriano se supone limitado por un nutriente clave. En contraste con un estudio anterior que mantuvo los niveles de neutrófilos y nutrientes constantes y requería términos de saturación para producir la biestabilidad, nuestro modelo muestra que los términos bilineales sencillos fomentan la biestabilidad cuando las densidades de nutrientes y neutrófilos pueden variar en respuesta a la densidad bacteriana. Un ejemplo aplicado a la bacteria Borrelia burgdorferi, que se alimenta de manganeso, ilustra por qué la respuesta rápida de los neutrófilos es clave para su capacidad de contener las infecciones bacterianas.

Palabras clave: biestabilidad; recursos limitados; fagocitosis; histéresis.

Mathematics Subject Classification: 92C37.

backward bifurcation in neutrophil-pathogen interaction 143

1 Introduction

Neutrophils are the most abundant type of white blood cell in most mammals. Highly mobile, with an average diameter of roughly 10 μm [5, 15], they can reach injured tissues faster than other types of cell. As part of the body’s nonspecific immune response, they can be sent to the site of an infection without the delays entailed in generating pathogen-specific lymphocytes, making them “first responders” against microbial infections. Neutrophils fight infection through phagocytosis, enveloping invading bacteria and decomposing them. Since the antimicrobial substances they produce during this process can also harm host tissues, macrophages then remove the neutrophils, which are short-lived in any case (mean lifetime less than a day in humans [22]). Neutrophil apoptosis also reduces inflammation at the infection site.

The nature of this immune response creates a complex relationship between neutrophils and the bacteria they fight, as the presence of the bacteria both stimulates their production and hastens their demise. In contrast, macrophages (which are larger and longer-lived) can attack many more pathogens before dying. Many models of immune response describe a predator-prey relationship between immune cells and pathogens. However, the dual nature of the influence of bacteria on neutrophil concentration sets it apart.

Previous studies have used mathematical models to describe interactions between microbial infections and nonspecific immune responses, often through the lens of a predator-prey relationship. Antia and Koella notably studied macrophage-parasite interaction in the presence of a constant population of complement factors (another element of the nonspecific immune system) [1]. With macrophages surviving phagocytosis longer to act on further particles, the macrophage concentration is affected only positively by the presence of the parasites, and the outcome of the interaction depends primarily on the phagocytosis rate, resulting in clearance, control, or merely stunting the growth of the pathogen. Similarly, Mayer et al. [14] presented a model of effector (immune) cells and target (pathogen) cells in which the effector cells are recruited in part by the presence of the pathogen (and in part by self-stimulation), but not destroyed by interacting with them. The rational polynomials used to model saturation in both immune cell recruitment terms lead to relatively complex behavior including bistability between pathogen extinction and survival, where pathogen survival can mean either an equilibrium concentration (control) or continued unbounded growth.

In 2000, Pilyugin and Antia incorporated handling time into the immune cell response to pathogens [17]. They developed a model in which the concentration of activated free immune cells X(t) is both increased (through recruitment) and decreased (through engagement in phagocytosis) by the presence of pathogens P [17, Equation 3.2], but in their analysis immediately passed to a simplified system [17, Equations 3.4 and 3.5] with a different form. While the original equation for X′(t) is linear in P, in the simplified system it is quadratic in P, with \frac{∂}{∂P} \frac{dX}{dt} negative for large enough P; that is, small infections have a net stimulatory effect on activated free phagocytes, but sufficiently large infections overwhelm the phagocytes by keeping them engaged while the remaining pathogen cells are free to replicate. The resulting analysis shows the bistability seen in [14] between either clearance or control, on the one hand, and continued unbounded growth of the pathogen on the other.

All of the above studies allow for unbounded pathogen growth by assuming no resource limitations. However, in many cases pathogen growth is limited by the availability of a key resource, often a metal such as iron or manganese. Iron is a significant mediator of host-pathogen interactions, as the body naturally sequesters it from pathogens and may reduce its availability even further when it detects infection [7]. Unbounded pathogen growth in a simple mathematical model can also violate the assumption of homogeneous mixing that underlies most descriptions of host-pathogen contact rates, necessitating a spatially explicit model.

One study which incorporated implicit resource limitations to pathogen growth was that of Smith et al. [21], which assumed logistic pathogen growth in studying innate immune responses to a pneumococcal lung infection. Their model included both neutrophils and multiple types of macrophages, with both positive and negative feedback from bacterial levels to neutrophil concentration. However, while the sacrificial nature of phagocytosis was modeled directly, the bacteria’s effect on neutrophil recruitment appears only indirectly: bacteria recruit alveolar macrophages, which increase cytokine production, which in turn recruit neutrophils to the infection site. The study’s analysis is purely numerical, but a bacterial survival threshold in the initial bacterial concentration suggests bistability, and the authors encourage further theoretical analyses in their conclusions.

Malka and colleagues [11, 12, 13] presented a model of bacterial growth in which both the bacteria’s resource (nutrient) and the neutrophil density were held constant, but saturation terms in both bacterial growth and bacteria-neutrophil interaction (attack rate) caused bistability.

The authors went so far as to assert that bacterial growth models without such saturation terms could never exhibit bistability [11].

On a larger scale, a qualitatively similar relationship to that between neutrophils and bacteria exists in predator-prey pairs where the prey uses a so-called lethal defense, typically via chemicals in its body that kill predators which eat it. For instance, a range of species from invertebrates such as gastropods, crabs, sea slugs, sea stars, ribbon worms, and blue-ringed octopi to pufferfish, roughskinned newts and certain types of frog generate the chemical tetrodotoxin [3]. This neurotoxin interferes with sodium channels in the brain, causing fatal paralysis. Similarly, monarch butterflies and some other butterfly species which feed on milkweed in the larval (caterpillar) stage sequester the cardenolides found in the plant [16], so that animals which eat the butterflies are then poisoned. However, studies of such populations typically focus on the evolutionary “arms race” in which predators develop resistance to each new prey defense mechanism.

This study proposes a more direct investigation of the dual effect of bacteria’s presence on neutrophil concentration, in an environment where bacterial growth is constrained by the availability of a key resource. The following sections develop and analyze a simple nonlinear dynamical system, offer a simple application to Borrelia burgdorferi using rates based in research literature, and conclude by interpreting the system’s behavior in the context of neutrophils’ role in fighting infections.

2 Model development

Neutrophil-pathogen interaction as described in the introduction involves a minimum of three interacting quantities: the concentrations of neutrophils N, bacteria B, and the limiting resource R on which bacterial growth depends. The model presented here, inspired by [4], consists of an ordinary differential equation to describe each quantity. The resource is assumed to be supplied at a constant rate θ from dietary and tissue sources, taken up by the host at a rate δ_{R}R, and used by the bacteria at a rate βBR proportional to the mass-action encounter rate at which bacteria find the resource. Bacterial growth is measured by multiplying the rate of resource uptake by a conversion efficiency parameter ϵ. Bacterial death due to encounters with neutrophils follows a mass-action encounter rate αBN, while bacterial death due to all other sources occurs at the linear rate δ_{B}B. This includes other immune system responses whose frequency (per capita rate) is independent of bacterial concentration.

Finally, neutrophils are assumed to be activated and recruited to the infection site at a baseline rate of ρ which is increased by γB when bacteria are present; they expire naturally at a rate δ_{N}N and undergo apoptosis following phagocytosis at a mass-action rate ωBN again proportional to the rate of neutrophil-bacteria encounters. This bilinear neutrophil death rate reflects how each phagocytosis event hastens the apoptosis of the neutrophil, building up bactericidal compounds which also signal macrophages for removal (but need not assume death after a single encounter). The resulting system is as follows:

\frac{dR}{dt} = θ − (βB + δ_{R})R,

\frac{dB}{dt} = (ϵβR − αN − δ_{B})B,  (1)

\frac{dN}{dt} = ρ + γB − (ωB + δ_{N})N,

with R measured in ng/mL, and B and N in cells/mL.

It should be noted that, in contrast to previous published models involving only a single variable, the bacterial concentration, this model uses negative feedback terms (not available to single-equation models) to incorporate limiting processes and saturation effects, rather than Monod-type functions. Simple bilinear (mass-action) terms involving two state variables produce the same limiting effects as saturation in an uncoupled equation does: that is, as bacterial concentration grows, it depletes the nutrient concentration, which in turn dampens bacterial growth. Similarly, as bacterial and neutrophil concentrations increase, their more frequent encounters reduce the other’s concentration, in turn reducing the encounter rate.

Parameter values used for this model are shown in Table 1 and are specific to the Lyme disease agent Borrelia burgdorferi. The lifespan of this bacteria ranges from several weeks to months [6]. Thus, assuming a lifespan of 4 weeks, we calculate the death rate of Borrelia burgdorferi to be 0.03 day^{−1}. Since Borrelia burgdorferi requires manganese (Mn) to grow [18], we take the limiting resource to be Mn. The estimated safe and adequate daily dietary intake (ESADDI) of Mn by adolescents and adults ranges from 2 to 5 mg [19].

However, only 1–5% of this is absorbed into the blood [2]. Since the human body has about 5000 mL of blood [8], we calculate the rate at which manganese is supplied to be at least 4 \frac{ng}{mL}. Assuming the normal concentration of Mn in the blood is 12 \frac{ng}{mL} (reference range 5–18 \frac{ng}{mL} [9]), back-calculation was employed (using the bacteria-free equilibrium described below) to arrive at a resource turnover rate of 0.33/day. The conversion efficiency parameter was adapted from the value presented in [4]. Since there is 16 times less Mn than iron (Fe) in BSK medium, (0.1 μM [23] of Mn compared to 1.6 μM of Fe [24]), we assume Borrelia burgdorferi is more efficient at utilizing Mn to replicate. The natural death rate of neutrophils, δ_{N}, is the reciprocal of the lifespan of neutrophils within the human body (taken from within the range given in [22]). Neutrophil production ρ is taken to be consistent with a baseline neutrophil concentration ρ/δ_{N} = 5.5 ∗ 10^{−6} cells/mL (taken from within the range given in [10]). Detailed information on the derivation of β, α, γ, ω can be found in [4].

Table 1: Model parameters and values.

Symbol Description  Value  Units Source
θ  Resource supply rate  4  \frac{ng}{mL∗day} [19]
δ_{R}  Resource turnover rate  0.33  \frac{1}{day}  this study
β  Bacteria growth affinity  2.85∗10^{−8} \frac{mL}{cell∗day} [4]
ϵ  Conversion efficiency parameter 1.6 ∗ 10^{8}  \frac{cells}{ng}  this study
α  Neutrophil attack rate  8.1 ∗ 10^{−5} \frac{mL}{cell∗day} [4]
δ_{B}  Bacteria death rate  0.03  \frac{1}{day}  [6]
ρ  Neutrophil recruitment rate  1.1 ∗ 10^{7}  \frac{cells}{mL∗day} [10]
γ  Neutrophil elicitation rate  0.8  \frac{1}{day}  [4]
δ_{N}  Neutrophil death rate  2  \frac{1}{day}  [22]
ω  Neutrophil apoptosis rate  1.8 ∗ 10^{−8} \frac{mL}{cell∗day} [4]

(following phagocytosis)

3 Analysis

In order to better understand the impact of neutrophils on the bacterial infection’s ability to persist, we preface our analysis of system (1) with a quick analysis of the subsystem without neutrophils, i.e., in which N = 0 (which then serves as a baseline for contrast). This two-dimensional system has a bacteria-free equilibrium R^{∗} = θ/δ_{R}, B^{∗} = 0, and an endemic equilibrium R^{∗} = δ_{B}/ϵβ, B^{∗} = \frac{δ_{R}}{β} (\frac{ϵβ}{δ_{B}} \frac{θ}{δ_{R}} − 1) which is positive iff R_{1} = \frac{ϵβ}{δ_{B}} \frac{θ}{δ_{R}} > 1.

The threshold quantity R_{1} is the bacterial reproductive number in the absence of neutrophils, the number of bacteria produced by each bacterium in the system before clearance. The expression can be interpreted as the ratio of the bacterial growth rate ϵβR^{∗} at the bacteria-free equilibrium, to the bacterial clearance rate δ_{B} in the absence of neutrophils. By linearization about each respective equilibrium, it can be seen that the bacteria-free equilibrium is locally asymptotically stable (LAS) iff R_{1} < 1, while the endemic equilibrium is LAS whenever it is positive (R_{1} > 1). The local stability in each case is seen to be global by applying Poincaré-Bendixson together with the observation^{1} that lim sup_{t→∞}(ϵR + B) ≤ ϵθ/ min(δ_{R}, δ_{B}), thus excluding unbounded solutions, and an application of Dulac’s Criterion with ϕ = 1/RB excluding periodic solutions in the positive quadrant. Thus R_{1} provides a complete measure of the bacterial infection’s ability to persist in the absence of neutrophils (it also figures in the analysis of the larger system, as will be seen next).

Returning to the full system (1), the bacteria-free equilibrium is R^{∗} = θ/δ_{R}, B^{∗} = 0, N^{∗} = ρ/δ_{N}. The endemic equilibrium condition can be written as a quadratic equation in B^{∗}, f(B^{∗}) = aB^{∗2} + bB^{∗} + c = 0, where

a = \frac{β}{δ_{R}} (δ_{B}ω + αγ) > 0,
b = \frac{β}{δ_{R}} (δ_{B}δ_{N} + αρ) + αγ + δ_{B}ω (1 − R_{1}) ,
c = δ_{B}δ_{N} + αρ − ϵβ \frac{θ}{δ_{R}} δ_{N} = (δ_{B}δ_{N} + αρ) (1 − R_{0}) ,
and R_{0} = \frac{ϵβ}{δ_{B} + α\frac{ρ}{δ_{N}}} \frac{θ}{δ_{R}} , R^{∗} = \frac{θ}{βB^{∗} + δ_{R}} , N^{∗} = \frac{ρ + γB^{∗}}{δ_{N} + ωB^{∗}} .

Here R_{0} is the bacterial reproductive number in the presence of neutrophils; by inspection R_{0} < R_{1}, with the denominator of R_{0} reflecting the additional bacterial removal rate due to phagocytosis produced by the baseline concentration (ρ/δ_{N}) of neutrophils.

Since a > 0, and R_{0} > 1 ⇔ c < 0, f has a unique positive root whenever R_{0} > 1. If instead R_{1} < 1 (which implies R_{0} < 1), all three coefficients are positive and f has no positive roots. In between, however, if R_{0} < 1 < R_{1}, then c > 0 but with (1 − R_{1}) < 0 in b, for ω large enough not only is b < 0 but (since a and b are linear in ω while c is independent of it) b^{2} − 4ac > 0, which implies that f has two positive roots. Since R_{0} is independent of ω, the conditions that R_{0} < 1 < R_{1} and ω be sufficiently large are independent.

1Since (ϵR + B)′ = ϵθ − δ_{R}ϵR − δ_{B}B ≤ ϵθ − (ϵR + B) min(δ_{R}, δ_{B}).

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 27(1): 141–156, Jan–Jun 2020

The specific condition on ω is^{2} ω > ω^{∗}, where

δ_{B}ω^{∗} = \frac{h + g}{R_{1} − 1} + 2\frac{h(1 − R_{0}) + \sqrt{h(1 − R_{0})[h(R_{1} − R_{0}) + gR_{1}(R_{1} − 1)]}}{(R_{1} − 1)^{2}}

and h = \frac{β}{δ_{R}} (δ_{B}δ_{N} + αρ), g = αγ.

It is straightforward to show using linearization that the bacteria-free equilibrium is LAS iff R_{0} < 1. If in addition R_{1} < 1, then one can show that the stability is in fact global, as follows. Suppose that B(0) > 0; then, for as long as B(t) > 0, \frac{dR}{dt} < θ − δ_{R}R, and thus R(t) < \frac{θ}{δ_{R}} + k_{1}e^{−δ_{R}t} for k_{1} = R(0) − \frac{θ}{δ_{R}}. Also R_{1} < 1 ⇔ ϵβ \frac{θ}{δ_{R}} < δ_{B}, so that

\frac{dB}{dt} < (ϵβ [\frac{θ}{δ_{R}} + k_{1}e^{−δ_{R}t}] − αN − δ_{B}) B < (ϵβk_{1}e^{−δ_{R}t} − αN) B. (2)

If k_{1} < 0, then dB/dt < 0, and B(t) → 0 as t → ∞. Otherwise, observe that \frac{dB}{dt} < ϵβk_{1}e^{−δ_{R}t}B, which implies B(t) < B(0) exp(k_{2} − e^{−δ_{R}t}) ≤ B(0) exp(k_{2}−1) for k_{2} = ϵβk_{1}/δ_{R}. Now dN/dt > ρ−(ωB+δ_{N})N > ρ−k_{3}N for k_{3} = ωB(0) exp(k_{2}−1)+δ_{N}, so that N(t) > \frac{ρ}{k_{3}} +(N(0) − \frac{ρ}{k_{3}} ) e^{−k_{3}t}. Ap-

plying this to (2) yields dB/dt<(−\frac{αρ}{k_{3}} + ϵβk_{1}e^{−δ_{R}t}−α (N(0) − \frac{ρ}{k_{3}}) e^{−k_{3}t}) B. From this last inequality it is clear that dB/dt < 0 beyond some time t^{∗}, and thus B(t) → 0 as t → ∞. This makes B a weak Lyapunov function for system (1), so the last step in proving global stability is to apply LaSalle’s invariant principle via the observation that on the set {(R, B, N) : B = 0} the bacteria-free equilibrium is the only attractor.

Although the stability analysis of any endemic equilibria is intractable here, numerical analysis indicates that the endemic equilibrium is asymptotically stable when unique, and when two exist, the second is unstable and generates a separatrix between the two basins of attraction.

2b = h + g + δ_{B}ω(1 − R_{1}) < 0 iff (δ_{B}ω) > (h + g)/(R_{1} − 1).

Then b^{2} − 4ac = (R_{1} − 1)^{2}(δ_{B}ω)^{2} − [2(h + g)(R_{1} − 1) + 4h(1 − R_{0})](δ_{B}ω)+ [(h − g)^{2} + 4R_{0}hg] > 0 iff (δ_{B}ω) is not between the roots of this latter quadratic, which are \frac{h+g}{R_{1}−1} +2^{\frac{h(1−R_{0})±\sqrt{h(1−R_{0})[h(R_{1}−R_{0})+gR_{1}(R_{1}−1)]}}{(R_{1}−1)^{2}}} . Since the radical in the numerator is greater than h(1 − R_{0}), the lower root is below (h + g)/(R_{1} − 1), so the only way to satisfy both b < 0 and b^{2} − 4ac > 0 is for (δ_{B}ω) to exceed the upper root.

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 27(1): 141–156, Jan–Jun 2020

Therefore, to summarize, system (1) has a unique endemic equilibrium iff R_{0} > 1, two endemic equilibria iff R_{0} < 1 < R_{1} and ω > ω^{∗}, and none otherwise. The bacteria-free equilibrium is LAS iff R_{0} < 1 and globally stable if R_{1} < 1. That is, through phagocytosis neutrophils reduce the bacterial load enough to eradicate sufficiently small initial infections (B(0) small enough that the initial condition falls within the bacteria-free equilibrium’s basin of attraction) which would otherwise persist (R_{0} < 1 < R_{1}), but if the sacrificial (apoptosis) aspect of phagocytosis is severe enough (ω sufficiently high), the neutrophils cannot control large initial infections (initial condition in an endemic equilibrium’s basin of attraction). This latter case corresponds to what is referred to as a backward bifurcation. Figure 1 superimposes the bifurcation diagrams of both system (1) and the corresponding (R, B) subsystem in a case where ω is high enough to cause a backward bifurcation. The shaded region in the first graph indicates the region of parameter space in which the presence of neutrophils leads to clearing the infection. For R_{0} in an interval below 1 (here, [0.5,1]), neutrophils are only able to clear sufficiently small initial infections.

(a) A stylized plot showing the region (shaded) of initial conditions for which the presence of neutrophils eradicates the infection. Without neutrophils, the infection dies out only to the left of the shaded region. To the right of the shaded region, the infection persists with or without neutrophils.

(b) A log plot showing the backward bifurcation for the B. burgdorferi–manganese parameter values in Table 1 but decreasing α by a factor of 5 and increasing ω by a factor of 50 (thus weakening the neutrophils), and using ϵ to vary R_{0}.

Figure 1: Superimposed bifurcation diagrams for system (1) and the corresponding (R, B) subsystem.

Substituting the parameter values for B. burgdorferi and manganese from Table 1 into the inequalities derived above reveals that in this scenario the bacteria are only prevented from invading by the presence of the neutrophils (R_{0} < 1 < R_{1}), and in addition the neutrophils are efficient and durable enough (ω < ω^{∗}) to preclude large endemic infections. However, if the neutrophils are weakened, the latter inequality can be reversed, leading to a backward bifurcation: if, for instance, α decreases by a factor of 5 and ω increases by a factor of 50, then ω > ω^{∗}, implying that although the neutrophils can stave off a small initial infection, they are unable to eradicate an established infection (see the second graph in Figure 1).

3.1 Comparison to Malka et al. [12]

The work of Malka and colleagues [11, 12, 13] presents a bacterial growth model consisting of a single equation,

B′(t) = ρ \frac{B}{1 + βB} + s − δ_{B}B − \frac{αNB}{1 + γB + ηN} ,  (3)

where the concentrations of neutrophils N and nutrients (implicit in ρ) are held constant, growth and neutrophil attack rate both saturate in the respective actors, and in addition a constant influx s of bacteria occurs (although in practice the authors set s to zero for analysis). In [11] the authors defined an entire class of similar models and wrote that models with growth and phagocytosis terms linear in B could not exhibit bistability. Our model and analysis show that replacing the saturation in these two terms with nutrient and neutrophil densities that respond to bacterial activity can lead to bistability. Both our model and equation (3) show that bistability requires terms representing interference with phagocytosis (γ in (3), ω in (1)) as well as resource limitations (β in both models). Both models also indicate the importance of the baseline neutrophil attack (phagocytosis) rate α in meeting bistability conditions.

Notably, analysis and discussion of (3) do not mention the bacterial infection’s basic reproduction number [11, 12, 13, 20], which is defined only when there is no constant influx, s = 0. In this case, the number is

R_{0} = \frac{ρ}{δ_{B} + α \frac{N}{1+ηN}} ,

whose structure closely parallels that for our model where η = 0 (ρ in (3) corresponds to ϵβθ/δ_{R} in (1), and N in (3) to ρ/δ_{N} in (1)).

In fact, the condition R_{0} < 1 for the expression above can be rewritten as

N > \frac{ρ}{α} (1 − \frac{δ}{ρ}) / (1 − \frac{(ρ − δ_{B})η}{α} ),

whose terms also appear in the condition the authors derive for bistability,

\frac{β}{γ} < (1 − \frac{δ}{ρ}) (1 − \frac{(ρ − δ_{B})η}{α} ).

The epidemiological perspective of our analysis describes the bistability directly in terms of reproduction numbers, and clarifies the role of neutrophils in this phenomenon.

4 Discussion

This study has shown how the “sacrificial” role that neutrophils play in fighting infection—a twist on the predator-prey relationship used to describe many immune cell–pathogen interactions—creates a complexity reflected in a simple mathematical model by a backward bifurcation and bistability. The relatively simple dynamical system used here illustrates how, if neutrophils are quick to arrive but of only mediocre efficacy (in terms of lifespan and pathogen encounter rate), they may be able to fight off a sufficiently small initial infection, but not a larger one; the simultaneous stability of both outcomes (eradication and endemic infection) also creates a hysteresis loop, making later eradication of a large enough initial infection much more difficult (as seen in the first graph in Figure 1, it may require cutting R_{0} in half). Nutrient limitations on bacterial growth, also play a key role in this bistability. Realistic parameter estimates for a B. burgdorferi infection dependent on manganese in the body, suggest an infection capacity in the right range for this to happen (R_{0} < 1 < R_{1}), but the neutrophils would have to be significantly weakened (perhaps in an immunocompromised patient) in order to meet the second and final criterion for backward bifurcation.

In contrast to previous studies demonstrating bistability in immune cell– pathogen interaction, this study illustrates how a simple (linear or bilinear), dynamic interaction among neutrophils, bacteria, and nutrients explains how neutrophils within a certain operating range can only clear sufficiently small bacterial infections. Those earlier studies assumed either (a) the possibility of unbounded bacterial growth [14, 17], a hypothesis which ignores pathogens’ dependence on limiting nutrients within a host, as well as the inaccuracy of a homogeneous mixing assumption for sufficiently large populations; or (b) constant levels of both nutrient and neutrophils [11, 12, 13].

The results of the present study are consistent with the behavior seen in Smith et al. [21] which suggested that bounded pathogen growth together with dual (positive and negative) feedback from pathogen to neutrophils may lead to a bacterial survival threshold including bistability. However, the model of Smith et al. also involves macrophages, and a roundabout feedback loop to neutrophils, whereas the present study uses a much simpler model with a direct feedback loop that makes clear the underlying mechanism. Our model also provides an alternative to the claim of Malka and colleagues, that growth and phagocytosis rates which do not explicitly saturate in bacterial concentration, cannot lead to bistability. The complex feedback between bacteria and neutrophil concentrations produces an environment which best eliminates infections when at low levels (unable to eradicate well-established infections), making neutrophils’ rapid response essential to their functionality.

To contextualize these results, it is worth noting that, in the human body, larger and longer-lasting infections are usually controlled by the body’s specific immune response, consisting of cells which typically survive numerous pathogen encounters. In contrast, neutrophils’ first-on-the-scene function is damage control until that specific immune response arrives, making the sacrificial nature of their encounters with pathogens (apoptosis and removal when they can envelop no more microbes) worth the reduced response time, as well as serving the function of reducing swelling at the infection site as time progresses.

The importance of nutrient limitations on bacterial growth increases in a context where different bacterial species infect the same host. Further work, already undertaken, considers how neutrophil-pathogen interaction impacts bacterial competition for such limited nutrients [4].

Acknowledgements

This project was inspired by [4], for which the authors thank their colleagues. CK dedicates this project to the memory of Fred and Carolyn Kribs.

References

[1] R. Antia, J.C. Koella, A model of non-specific immunity, Journal of Theoretical Biology 168(1994), no. 2, 141–150. doi: 10.1006/ jtbi. 1994.

    1094

[2] M.G. Baker, C.D. Simpson, B. Stover, L. Sheppard, H. Checkoway, B.A. Racette, N.S. Seixas, Blood manganese as an exposure biomarker: state of the evidence, Journal of Occupational and Environmental Hygiene

    11(2014), no. 4, 210–217. doi: 10.1080/ 15459624. 2013. 852280

[3] V. Bane, M. Lehane, M. Dikshit, A. O’Riordan, A. Furey, Tetrodotoxin: chemistry, toxicity, source, distribution and detection, Toxins 6(2014),

no. 2, 693–755. doi: 10.3390/ toxins6020693

[4] A. Carter, A. Toja, O. Olawoyin, J. Gonzalez, J. Grover, H. Kojouharov, C. Kribs Zaleta, A theoretical model of coinfection dynamics: modeling competition dynamics between Borrelia burgdorferi and Anaplasma phagocytophilum within a human host, Mathematics Department, Preprint Series 2014-06, University of Texas at Arlington, 2014.

[5] V.G. Daniels, P.R. Wheater, H.G. Burkitt, Functional histology: A text and colour atlas, Churchill Livingstone, Edinburgh, 1979.

[6] A.N. Datar, N. Kaur, S. Patel, D.F. Luecke, E. Sapi, In vitro effectiveness of samento and banderol herbal extracts on the different morphological forms of Borrelia burgdorferi, Ph.D. Thesis, Univ. of New Haven, West Haven CT, 2010.

[7] H. Drakesmith, A.M. Prentice, Hepcidin and the iron-infection axis, Science 338(2012), no. 6108, 768–772. doi: 10.1126/ science.

1224577

[8] L. Geggel, How much blood is in the human body, Live-

Science, 2016. Available at https:// www. livescience.  com/

32213-how-much-blood-is-in-the-human-body.      html

[9] R.B. Jain, Y.S. Choi, Normal reference ranges for and variability in the levels of blood manganese and selenium by gender, age, and race/ethnicity for general U.S. population, Journal of Trace Elements in Medicine and Biology 30(2015), 142–152. doi: 10.1016/ j.jtemb. 2014. 12.004

[10] Y. Li, A. Karlin, J.D. Loike, S.C. Silverstein, A critical concentration of neutrophils is required for effective bacterial killing in suspension, Proceedings of the National Academy of Sciences of the United States of America 99(2002), no. 12, 8289–8294. doi: 10.1073/ pnas.

122244799

[11] R. Malka, V. Rom-Kedar, Bacteria-phagocyte dynamics, axiomatic modelling and mass-action kinetics, Mathematical Biosciences and Engineering 8(2011), no. 2, 475–502. doi: 10.3934/ mbe. 2011. 8.475

[12] R. Malka, E. Shochat, V. Rom-Kedar, Bistability and bacterial infections, PLoS ONE 5(2010), no. 5, 1–10. doi: 10.1371/ journal. pone.

0010010

[13] R. Malka, B. Wolach, R. Gavrieli, E. Shochat, V. Rom-Kedar, Evidence for bistable bacteria-neutrophil interaction and its clinical implications, Journal of Clinical Investigation 122(2012), no. 8, 3002–3011. doi: 10.

1172/ JCI59832

[14] H. Mayer, K. S. Zaenker, U. an der Heiden, A basic mathematical model of the immune response, Chaos: An Interdisciplinary Journal of Nonlinear Science 5(1995), no. 1, 155–161. doi: 10.1063/ 1.166098

[15] M.J. Niemiec, B. De Samber, J. Garrevoet, E. Vergucht, B. Vekemans, R. De Rycke, E. Björn, L. Sandblad, G. Wellenreuther, G. Falkenberg, P. Cloetens, L. Vincze, C.F. Urban, Trace element landscape of resting and activated human neutrophils on sub-micrometer level, Metallomics 7(2015), no. 6, 996–1010. doi: 10.1039/ c4mt00346b

[16] G. Petschenka, A. A. Agrawal, Milkweed butterfly resistance to plant toxins is linked to sequestration, not coping with a toxic diet, Proceedings of the Royal Society B 282(2015), no. 1818: 20151865, 1–9. doi: 10.1098/

rspb. 2015. 1865

[17] S.S. Pilyugin, R. Antia, Modeling immune responses with handling time, Bulletin of Mathematical Biology 62(2000), no. 5, 869–890. doi: 10.

1006/ bulm. 2000. 0181

[18] J.E. Posey, F.C. Gherardini, Lack of a role for iron in the Lyme disease pathogen, Science 288(2000), no. 5471, 1651–1653. doi: 10.1126/

science. 288. 5471. 1651

[19] A.B. Santamaria, Manganese exposure, essentiality & toxicity, Indian Journal of Medical Research 128(2008), no. 4, 484–500.

[20] S.C. Silverstein, R. Rabadan, How many neutrophils are enough (redux, redux)?, Journal of Clinical Investigation 122(2012), no. 8, 2776–2779.

doi: 10.1172/ JCI63939

[21] A M. Smith, J.A. McCullers, F.R. Adler, Mathematical model of a threestage innate immune response to a pneumococcal lung infection, Journal of Theoretical Biology 276(2011), no. 1, 106–116. doi: 10.1016/ j.

jtbi. 2011. 01.052

[22] T. Tak, K. Tesselaar, J. Pillay, J.A. Borghans, L. Koenderman, What’s your age again? Determination of human neutrophil half-lives revisited, Journal of Leukocyte Biology 94(2013), no. 4, 595–601. doi: 10.1189/ jlb.

1112571

[23] B. Troxell, M. Ye, Y. Yang, S.E. Carrasco, Y. Lou, X. F. Yang, Manganese and zinc regulate virulence determinants in Borrelia burgdorferi, Infection and Immunity 81(2013), no. 8, 2743–2752. doi: 10.1128/ IAI.

00507-13

[24] J.R. Zhang, S.J. Norris, Kinetics and in vivo induction of genetic variation of vlsE in Borrelia burgdorferi, Infection and Immunity 66(1998), no. 8, 3689–3697.