Bistable Mathematical Model of Neutrophil Migratory Patterns After LPS-Induced Epigenetic Reprogramming

The highly controlled migration of neutrophils toward the site of an infection can be altered when they are trained with lipopolysaccharides (LPS), with high dose LPS enhancing neutrophil migratory pattern toward the bacterial derived source signal and super-low dose LPS inducing either migration toward an intermediary signal or dysregulation and oscillatory movement. Empirical studies that use microfluidic chemotaxis-chip devices with two opposing chemoattractants showed differential neutrophil migration after challenge with different LPS doses. The epigenetic alterations responsible for changes in neutrophil migratory behavior are unknown. We developed two mathematical models that evaluate the mechanistic interactions responsible for neutrophil migratory decision-making when exposed to competing chemoattractants and challenged with LPS. The first model, which considers the interactions between the receptor densities of two competing chemoattractants, their kinases, and LPS, displayed bistability between high and low ratios of primary to intermediary chemoattractant receptor densities. In particular, at equilibrium, we observe equal receptor densities for low LPS (< 15ng/mL); and dominance of receptors for the primary chemoattractant for high LPS (> 15ng/mL). The second model, which included additional interactions with an extracellular signal-regulated kinase in both phosphorylated and non-phosphorylated forms, has an additional dynamic outcome, oscillatory dynamics for both receptors, as seen in the data. In particular, it found equal receptor densities in the absence of oscillation for super-low and high LPS challenge (< 0.4 and 1.1 376 ng/mL). Predicting the mechanisms and the type of external LPS challenge responsible for neutrophils migration toward pro-inflammatory chemoattractants, migration toward pro-tolerant chemoattractants, or oscillatory movement is necessary knowledge in designing interventions against immune diseases, such as sepsis.


INTRODUCTION
Researchers have recently challenged the dogma that innate immunity is the same at every challenge. It has been shown that macrophages are able to develop different kinds of memory depending on the type of priming they encounter via epigenetic reprogramming (Yuan et al., 2016a,b). For instance, they can develop a memory phenotype that leads them to be less reactive or even tolerant to a challenge, or they can develop a memory phenotype that leads them to have an enhanced response to a challenge. This same concept has recently been shown by us for neutrophil migratory decision-making, and it is thought that the response is influencing the outcomes of infectious diseases. For example, in sepsis or COVID-19 infection, the immune system overreacts because of underlying low-grade inflammation that primes neutrophils into choosing between tolerant and inflammatory migratory phenotypes (Alves-Filho et al., 2005. As a result, neutrophils can migrate to healthy organs and unleash their anti-microbial arsenal in healthy tissue, leading to organ failure in the lungs, kidney, or heart. The mechanisms underlying trained innate immunity have not been fully elucidated, with epigenetic modifications playing a key role in the induction of innate memory or training (Pillay et al., 2010;Demaret et al., 2015). In this study we investigate innate memory in the context of neutropil migratory decision-making.
The ability of neutrophils to migrate plays a pivotal role in a cell's ability to clear infections and resolve inflammation. During infection and inflammation, chemoattractants are released, signaling and activating neutrophils in the bloodstream. Neutrophils must be able to precisely migrate within the tissue to the specific site of infection, without being diverted toward other locations, in a process called chemotaxis. Chemotaxis is a highly regulated process that involves activation of various pathways and downstream polarization of the cell (Kolaczkowska and Kubes, 2013). The first step in chemotaxis is recognition of chemoattractants by the cell. Cells have specific receptors on their surface for various chemoattractants. These chemoattractant receptors are G protein-coupled receptors (GPCRs), which are regulated by a variety of G protein-coupled receptor kinases (GRKs) (Murphy, 1994;Dianqing, 2005). When bound by a specific agonist, in this case a chemoattractant, the GPCRs undergo phosphorylation, which unbinds the G proteins and desensitizes the receptor. This leads to internalization of the receptor, activation of downstream signaling pathways, and activation of cellular responses, such as cell polarization and chemotaxis (Murphy, 1994;Dianqing, 2005;Futosi et al., 2013). After internalization, receptors can be recycled back to the cell surface, where they can again be bound by the receptor's agonist. This process is crucial in chemotaxis, as it allows the cell to continue sensing the chemoattractant and migrate in its direction (Neel et al., 2005). Most chemoattractant receptors are similar in their response to ligand-binding; however there are slight differences in the activated signaling pathways (Heit et al., 2002(Heit et al., , 2008. Within the tissue, neutrophils are exposed to several chemoattractants at once, originating from pathogens, cells within the tissue, the endothelium, and several other sources (Kolaczkowska and Kubes, 2013). Cells must prioritize these signals to properly clear the pathogen. It has been hypothesized that neutrophils have an internal hierarchy, where chemoattractants derived from bacterial sources and the complement system, such as fMLP and C5a (Heit et al., 2002;Petri and Sanz, 2018), take precedent over intermediary chemoattractants, such as LTB 4 and IL-8, which are secreted by other immune cells. This leads to neutrophils migrating toward end-target chemoattractants over intermediary chemottractants in a competitive environment (Heit et al., 2002(Heit et al., , 2008Wang et al., 2016b), allowing neutrophils to prioritize an invading pathogen. This hierarchy is thought to occur through the activation of differing signaling pathways, where end-target chemoattractants signal through p38 MAPK and intermediary chemoattractants signal through PI3K (Heit et al., 2002(Heit et al., , 2008. The highly controlled migration of neutrophils toward the site of an infection, as well as their dynamic interaction with pathogens, can be altered when they are pre-conditioned with Lipopolysaccharides (LPS) to induce endotoxin priming. In previous work, we showed that training with high dose LPS (100 ng/mL) enhances neutrophil migration toward the endtarget, bacterial derived, source signal fMLP. By contrast, training with super-low dose LPS (1 ng/mL) alters neutrophil migratory phenotypes, which either migrate toward the intermediary signal LTB 4 or become dysregulated and exhibit oscillatory migratory patterns (Jones et al., 2016;Boribong et al., 2019). While the empirical data shows that neutrophils trained with LPS change migratory phenotype, it does not give information on the molecular mechanisms responsible for the difference in behavior. The migratory decision-making process is finely governed by complex signaling networks that dynamically receive and interpret molecular and cellular signals from outside and within. The intrinsic complexity of immune cell decision-making processes has created difficulty for experimental immunologists to determine the mechanisms of disease, in spite of expansive experimental studies with conventional reductionist cellular and molecular approaches. It is increasingly recognized that crossdisciplinary studies combining experimental and mathematical modeling approaches are critically required.
In this study, we investigate the molecular mechanisms of neutrophil migratory decision-making in the presence of competing chemoattractants and external challenge with LPS, by building deterministic mathematical models of interaction between two chemoattractant receptors, Formyl Peptide Receptor 1 (FPR1) and Leukotriene B4 Receptor 1 (BLT1), and key molecules involved in their regulation. We are interested in determining the relationship between the receptor dynamics and migration pattern, and in quantifying the LPS dose resulting in neutrophils migration toward a pro-inflammatory chemoattractant, toward a pro-resolution chemoattractant, or in neutrophils dysregulation and oscillation (Fan and Malik, 2003;Liu et al., 2012;Byrne et al., 2014). The model will qualitatively match the experimental results of our previous work, where stimulation with a super-low concentration of LPS will result in greater BLT1 over FPR1, and stimulation with a high concentration of LPS will result in greater FPR1 over BLT1 . We construct a model with bistable behavior, with the motif for bistability coming from the non-linear mutual inhibition of GRK2 and GRK5 (see Figure 1). The dual inhibition leads to the activation of different signaling pathways (p38/JNK vs. ERK), leading to differences in functional neutrophil migration (Davenport et al., 2020). Both GRK2 and GRK5 have been demonstrated to be critical mediators of the molecular alterations that occur in the inflammatory disorders, but the complex mutual inhibition interaction has largely been ignored (Philipp et al., 2014). Mathematical models have been used before to model cellular decision-making Kadelka et al., 2019), neutrophil chemotaxis (Ionides et al., 2004;Postma and van Haastert, 2016;Bayani et al., 2020), immune responses Fischer, 2008;Nelson et al., 2009;Vodovotz et al., 2009) and bistable dynamics (Ciupe et al., 2007(Ciupe et al., , 2018Leber et al., 2016).

Mathematical Model of Migratory Decision-Making
We developed a novel system of differential equations based on diagram in Figure 2 (Prossnitz et al., 1995;Arraes et al., 2006;Sorriento et al., 2008;Wang et al., 2016a). For simplicity, we model linear effects of LPS on the kinases' activity. In particular, we assume that the GRK2 activation occurs at rate c w + a w [LPS], with c w and a w being the LPS-independent and LPS-dependent activation rates. Similarly, GRK5 activation occurs at rate c f + a f [LPS], with c f and a f being the LPS-independent and LPS-dependent activation rates. The two kinases mutually inhibit one another. We model inhibition of GRK2 via GRK5 at rate 1/(b fw + [GRK5] n ) and inhibition of GRK5 via GRK2 at rate 1/(b wf + [GRK2]), where b fw and b wf are the mutual inhibition rates of GRK2 by GRK5 and GRK5 by GRK1, respectively. n is the cooperativity coefficient. We assumed increased cooperativity in GRK2 inhibition by GRK5, but not the inhibition of GRK5 by GRK2. The results are preserved if the same cooperativity is included in the GRK5 inhibition by GRK2 (not shown). We assume GRK2 and GRK5 decay at per capita rates d w and d f , respectively, with GRK5 decay being modeled in a density dependent manner, with the GRK5 value where the decay is half-maximal being given by parameter b f .
The chemoattractant receptors FPR1 and BLT1 internalize from the plasma membrane into the cell via phosphorylation (Magalhaes et al., 2012;Mócsai et al., 2015). We assume that the number of receptors on a cell is conserved and, through the process of dephosphorylation, the receptors are recycled and brought back to the surface of the cell. Thus, we have conservation laws of the total number of the receptor equalling the sum of the non-phosphrylated and phosphorylated receptor, The process of receptor phosphorylation and dephosphorylation is modeled using Hill-type functions. In particular, FPR1 is produced through dephosphorylation, modeled by a Michaelis-Menten where a 1 is maximal production and J F 1 is the receptor quantity where dephosphorylation is half-maximal. Similarly, FPR1 is lost through phosphorylation, which is enhanced in the presence of GRK2 (Wang et al., 2016a). We model this by a Hill-type function where a 2 is the maximal rate and J F 2 is the receptor quantity where phosphorylation is half-maximal.
BLT1 is produced through dephosphorylation, modeled by a Michaelis- where b 1 is the maximal production rate and J B 1 is the receptor quantity where dephosphorylation is half-maximal. BLT1 is lost through phosphorylation, which is enhanced in the presence of both GRK2 and GRK5 (Gaudreau et al., 2002;Chen et al., 2004). We model this by a Hill-type where b 2 are b 3 are maximal decay rates and J B 2 is the receptor quantity where phosphorylation is half-maximal. We assume a single LPS dose, after which LPS decays exponentially at a rate d L (Kadelka et al., 2019). The dynamical system describing these interactions is given by: (1) We are interested in determining the ratio between the cells that migrate toward the primary and those that migrate toward the intermediary chemoattractants given, as a proxy, by the ratio of their receptors FPR1/BLT1, when initial LPS is varied.

Experimental Data
In previous research, we used a microfluidic competitive chemotaxis-chip device to measure the migratory decisionmaking process of dHL-60 cells, a model neutrophil cell line, 5 h after they were pre-challenged with super-lowdose (1 ng/mL) and high-dose (100 ng/mL) of LPS in the presence of two competing chemoattractants, LTB4 and fMLP . Challenging the cells with a super-low dose of LPS resulted in fMLP/LTB4 ratio of 0.8672. Challenging the cells with a high dose of LPS (100 ng/mL) resulted in fMLP/LTB4 ratio of 10.2646 (see Figure 3A).

Experimental Data
Experimental results reported that neutrophils treated overnight with LPS may lose their ability to move up the chemoattractant gradient, become disoriented, and display oscillatory behavior . Moreover, the highest number of cells to display such oscillatory behavior occurs following LPS exposure with super-low dose (1 ng/mL) (see Figure 3B) .

Bistable FPR1 and BLT1 Dynamics
We evaluated neutrophil migration between end-target chemoattractant fMLP and intermediary chemoattractant LTB 4 by developing model (1) Figure 4B). Five hours following challenge with super-low-dose (1 ng/mL) LPS, model (1) predicts the presence of a small number of receptors, which are distributed equally among FPR1 and BLT1, [FPR1]/[BLT1](5) = 1 (see Figure 4A, solid lines). Under our abstraction this means that, following challenge with 1 ng/mL LPS, an equal number of neutrophils migrated toward the fMLP and LTB4 gradients. Conversely, 5 h following highdose challenge (100 ng/mL) LPS, model (1)  To determine the relationship between the LPS challenge dose and the FPR1/BLT1 ratio, we derived a graph that quantifies [FPR1]/[BLT1](5), 5 h following cell priming, as a function of the [LPS] dose, predicted by model (1) and parameter values/initial conditions in Table 1. We found that the experimental observation for the super-low-dose (1 ng/mL) LPS, [FPR1]/[BLT1](5) = 1, is preserved for all challenges with LPS values lower than 3.9 ng/mL. For 4 − 6.7 ng/mL LPS challenge, [BLT1] exceeds [FPR1] at t = 5 h, but the two receptors will eventually reach identical levels at equilibrium.  (1). Parameters and initial conditions are given in Table 1. Lastly, the [FPR1]/[BLT1](5) ratio grows larger than one and keeps increasing for LPS dose > 6.7 ng/mL, eventually reaching the experimental prediction of ten, [FPR1]/[BLT1](5) = 10, for high-dose LPS challenge (100 ng/mL) (see Figure 5) and increasing further as time passes or for higher challenge (not shown).

Long-Term Results and Motifs of Bistability
We have chosen the parameters in model (1)

Molecular Mechanisms of Cell Oscillatory Migration
Experimental results reported that neutrophils treated overnight with LPS may loose their ability to move up the chemoattractant gradient, become disoriented, and display oscillatory behavior . Moreover, the highest number of cells to display such oscillatory behavior occurs following LPS exposure with low dose of 1 mg/ml (see Figure 3B) . To determine the molecular mechanisms responsible for the oscillations, we extended the bistable system (1), by coupling it with an activator-inhibitor oscillatory model for the dynamics of non-phosphorylated and phosphorylated extracellular signalregulated kinases, [ERK] and [ERK p ], and two auxiliary enzymes, [E] and [E p ] based on diagram (7), model (4) and parameter values/initial conditions in Tables 1, 2. Moreover LPS is fixed at its initial condition. Under the chosen parameters, we obtain long-term oscillatory movement for all populations, for constant super-low-challenge LPS (1 ng/mL), as predicted by the data ) (see Figures 8A-C). Interestingly, the oscillatory behavior is maintained for a short LPS range, (0.5-1.  (1), the oscillatory dynamics are due to the oscillatory dynamics of the [ERK] and [ERKp] kinases. This oscillatory behavior can be broken by either increasing or decreasing the constant LPS challenge. Such information can inform interventions, as dysoriented neutrophil movement is not desirable and has been shown to have negative effects during pathogenic infections. Dysregulated neutrophil response to infection can lead to sepsis and end-organ failure and is a leading cause of death worldwide (Reddy and Standiford, 2010;Shen et al., 2017).

DISCUSSION
In this study, we developed compartmental mathematical models of molecular interactions that govern neutrophil migratory FIGURE 7 | Diagram for model (4). patterns when exposed to competing chemoattractants and challenged with external stimuli. When the models were restricted to the interactions between the chemoattractants' receptors, their kinases, and LPS, we predicted a bistable switch between two states: one in which the densities of the two chemoattractant receptors, FPR1 and BLT1, are equal and one in which the receptors for the primary chemoattractant, FPR1, dominate. We hypothesized that the two states correspond to two states observed experimentally: equal migration toward the primary and intermediary chemoattractants, fMLP and LTB4, and predominant migration toward the primary chemoattractant, fMLP . The experimental data connected the differential migratory outcomes with the magnitude of the external LPS challenge, with super-low 1 ng/mL LPS leading to equal migration toward both chemoattractants and high 100 ng/mL LPS leading to ten times higher migration toward fMLP, 5 h after challenge. In the mathematical model, the external signal corresponds to the initial condition for variable LPS. Our model was calibrated to match the experimental data 5 h after stimuli, with [FPR1]/[BLT1](5) = 1 for [LPS](0) = 1 and [FPR1]/[BLT1](5) = 10 for [LPS](0) = 100. Furthermore, 5 h after LPS challenge, we obtain equal levels of FPR1 and BLT1 receptors for all initial conditions [LPS](0) < 3.7 ng/mL and increasingly more FPR1 than BLT1 receptors when the initial condition for LPS increases, with ten times more FPR1 than BLT1 receptors when [LPS](0) = 100 ng/mL. When we run the model to equilibrium, however, we obtain equal FPR1 and BLT1 receptors for an even larger range of LPS initial conditions, < 15 ng/ml; and dominant FPR1 levels for LPS> 16 ng/mL. This implies that not just the challenge dose, but the duration of the experiment may influence the quantitative outcomes.
We investigated the molecular interactions that are responsible for the bistable outcomes and found that when the mutual inhibition of GRK2 and GRK5 kinases is either removed, or the balance is broken, the model is no longer bistable. Instead, when run to equilibrium, it settles into a state where FPR1 receptors dominate the outcome, indicative of predominant migration toward the primary chemoattractant.
When the models were expanded to add the interaction with the ERK signaling pathways under constant LPS challenge, we obtained a third dynamical state, where we have equal FPR1/BLT1=1, but the equilibrium is lost, and the FPR1 and BLT1 receptors oscillate between two values. We assume this to be indicative of neutrophil oscillation which, in the experimental setup, is equivalent to cells changing direction . We investigated how changes in the constant LPS dose affect the outcomes and found that, at equilibrium, FPR1 and BLT1 receptors are equal and non-oscillating for both supersuper-low (< 0.5 ng/mL) and high dose (1.1, 375 ng/mL) LPS; are equal and oscillating for super-low dose (0.5-1.1 ng/mL) LPS; and FPR1 outnumbers BLT1 for super-high dose (> 376 ng/mL) LPS. The non-asymptotic dynamics are due to the oscillatory behavior of phosphorylated and non-phosphorylated ERK molecules, who undergo an auto-catalytic interaction with two undefined enzymes. We have modeled the interaction using an activation-inhibition motif, with similar dynamics being  obtained if the oscillations are induced by a substrate depletion motif (Tyson et al., 2003). Further work is needed to determine the nature of enzyme or their regulation. We are currently working to validate the model by retrieving neutrophils from our microfluidic device post-migration and quantifying FPR1, BLT1, GRK2, and GRK 5 levels by Droplet Digital TM PCR (ddPCR TM ).
Our models are limited by the presence of many unknown parameters. While we strived to match the empirical data, the results are mostly qualitative. Similar results can be obtained with many different parameters sets. For example, c w only slightly influences the FPR1/BLT1 ratio at t = 5 h (see Figure 9, left panel), while the cooperativity coefficient n and the LPSdependent GRK5 production a f have drastic effects on the size of the ration (see Figure 9, middle and right panels). The unchanging factors, however, are the motifs of bistability, which are induced by the dual inhibition of the G-protein kinases; the oscillatory motifs, which are induced by the oscillatory ERK dynamics; and the influence of external stimuli on outcomes.
In conclusion, we developed mathematical models for the molecular interactions responsible for neutrophils migratory phenotypes, calibrated them against empirical data, and used their dynamics to determine the external stimuli ranges that account for neutrophils migration toward a pro-inflammatory chemoattractant, a pro-tolerant chemoattractant, or oscillatory dynamics indicative of dysorientation and loss of function. Understanding the relationship between neutrophils' dynamics and the mechanisms responsible for their movement is important for preventing and predicting immune disorders. Activation markers in neutrophils are potential biomarkers for the diagnosis and prognosis of sepsis. Septic patients can be screened for neutrophil markers, such as GRK 2/5, FPR1, and BLT1, and this predictive model can guide patient treatment. In the future, we plan to expand this model to include more complex signaling from the microenvironment and aim to predict not only dysfunctional migration in neutrophils, but even the probability of cells accumulating in specific organs, such as the lung or kidney. We can use these predictive models to define optimal patient treatment and to identify immunotherapeutic targets (i.e., small molecule inhibition, microRNAs, gene therapy) to promote directional neutrophil migration.

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