A clinically parameterized mathematical model of Shigella immunity to inform vaccine design

We refine and clinically parameterize a mathematical model of the humoral immune response against Shigella, a diarrheal bacteria that infects 80-165 million people and kills an estimated 600,000 people worldwide each year. Using Latin hypercube sampling and Monte Carlo simulations for parameter estimation, we fit our model to human immune data from two Shigella EcSf2a-2 vaccine trials and a rechallenge study in which antibody and B-cell responses against Shigella′s lipopolysaccharide (LPS) and O-membrane proteins (OMP) were recorded. The clinically grounded model is used to mathematically investigate which key immune mechanisms and bacterial targets confer immunity against Shigella and to predict which humoral immune components should be elicited to create a protective vaccine against Shigella. The model offers insight into why the EcSf2a-2 vaccine had low efficacy and demonstrates that at a group level a humoral immune response induced by EcSf2a-2 vaccine or wild-type challenge against Shigella′s LPS or OMP does not appear sufficient for protection. That is, the model predicts an uncontrolled infection of gut epithelial cells that is present across all best-fit model parameterizations when fit to EcSf2a-2 vaccine or wild-type challenge data. Using sensitivity analysis, we explore which model parameter values must be altered to prevent the destructive epithelial invasion by Shigella bacteria and identify four key parameter groups as potential vaccine targets or immune correlates: 1) the rate that Shigella migrates into the lamina propria or epithelium, 2) the rate that memory B cells (BM) differentiate into antibody-secreting cells (ASC), 3) the rate at which antibodies are produced by activated ASC, and 4) the Shigella-specific BM carrying capacity. This paper underscores the need for a multifaceted approach in ongoing efforts to design an effective Shigella vaccine.

Introduction rather than individual level. We use Latin hypercube sampling, least squares optimization, and Monte Carlo simulations to explore the relevant multi-dimensional parameter space and to identify best-fit parameter value combinations. We also update our mathematical model introduced in [29] with key structural changes. Most notably, we use ordinary instead of delay differential equations for numerical efficiency while carrying out Monte Carlo runs during parameter estimation. We also directly include memory B cell stem-like activity (that is, intermittent asymmetric division of B M to maintain the B M pool) to counteract the degree of memory cell depletion as B M differentiate into ASC that was observed with the first model, which assumed a baseline logistic B M growth rate and then only purely symmetric differentiation of B M into ASC. After comprehensive clinical parameterization, we perform identifiability analysis and use numerical sensitivity analysis to identify the few critical parameter values that govern the clearance and severity of Shigella infections. With these, we make model-based predictions about which biological interactions might-and which biological interactions likely do not-correlate with protective immunity against Shigella. We find that antibody-based vaccines targeting Shigella 0 s outer membrane components such as LPS and OMP are unlikely to elicit sufficient immunity for protection at a group level except possibly with major boosts in particular immune parameters, which we identify. Our modeling suggests secondary humoral immune targets upon which to focus in future vaccine development efforts. These include 1) the rate that Shigella migrates into the lamina propria or epithelium, 2) the rate that B M differentiate into ASC, 3) the rate at which antibodies are produced by activated ASC, and 4) the Shigella-specific B M carrying capacity (that is, the number of Shigella-specific B M at homeostasis between infections). Future modeling and vaccine efforts could also examine the efficacy of CMI and other immune targets not directly measured in these clinical studies and thus not included in the present model.

Data
We use clinical data supplied and published by our collaborators from three human clinical studies in the 1990s: two vaccine clinical trials followed by a rechallenge study [30][31][32][33][34]. The two vaccine trials (called the EcSf2a-2 studies in this paper or ET in labels) tested the efficacy of EcSf2a-2, a hybrid E. coli-Shigella flexneri 2a vaccine candidate [30][31][32]. This clinical trial separated volunteers into control and vaccinated groups, the latter of which were given the EcSf2a-2 vaccine. All volunteers were later challenged with 1,400 colony forming units (cfu) of 2457T wild-type (wt) Shigella. Volunteers were tracked for a month and time-course data were gathered at days 0, 7, and 10 for ASC of IgA-and IgG-type and at days 0, 14, 21, and 28 for antibodies-specifically, serum IgA and serum IgG directed at LPS and OMP in Shigella 0 s outer membrane. In the EcSf2a-2 vaccine studies, 14 volunteers received placebo while 45 volunteers received the vaccine in three or four doses. All volunteers then received the wild-type challenge with 2457T Shigella. We aggregate the respective volunteers in the EcSf2a-2 trials to form the placebo+wt primary infection data, which we label ET: 2457Tx1 for the EcSf2a-2 studies with no vaccine and one challenge of 2457T Shigella, and separately the vaccine+wt secondary infection data, which we label ET: EcSf2a-2x1, 2457Tx1 for the EcSf2a-2 trials with one round of EcSf2a-2 vaccine and then one challenge of 2457T vaccine. (One vaccine round includes all three or four doses, as applicable.) We fit our model against the mean ASC and antibody levels for each group.   [30][31][32]34]. Measured antibody (serum IgA, serum IgG), IgA-secreting ASC (A-ASC), and IgG-secreting ASC (G-ASC) directed against lipopolysaccharide (LPS, left) and O-membrane proteins (OMP, right) are given over time for each trial. We convert measured specific antibody titers to the number of specific Ab/mL by a factor of 1: 1.204 × 10 10 Ab/mL, as described in the text. Measurements were taken 0, 10, 21, and 28 days after challenge for antibodies and on days 0, 7, and 10 for ASC; displays each day are In a follow-up to the EcSf2a-2 trials, 30 volunteers were challenged or rechallenged with 1,400 cfu of 2457T wt Shigella. We call this the 2457T study or 2T in labels. The 19 formerly immunologically naive individuals who received a 2457T challenge are labeled 2T: 2457Tx1, for the 2T study with one challenge of wt 2457T Shigella. Six other individuals, labeled 2T: 2457Tx2, had served as controls in the EcSf2a-2 study; these individuals had previously been challenged with 2457T Shigella without receiving the EcSf2a-2 vaccine candidate and in the second study were challenged a second time with wt 2457T Shigella. These data sets, also shown in Fig 1, thus show antibody and ASC levels directed against LPS and OMP both for formerly naive individuals after a primary infection and for previously exposed individuals after a secondary wild-type infection. In theory, the latter individuals should exhibit the reactivation of any natural immunity that was established during the primary infection against Shigella, although the level of immunity first established was not necessarily protective. A third category of 5 volunteers who received the EcSf2a-2 vaccine and challenge in the first trial followed by rechallenge in the second study (thus, possibly exhibiting tertiary immunity) are labeled 2T: EcSf2a-2x1, 2457Tx2 and are not included in this paper's analysis to avoid mixing the effects of vaccine-induced versus wt-induced immunity.
Notice in Fig 1 that from primary to secondary immune responses there is in many cases no statistical rise and often slight declines in the levels of serum Ab and ASC directed against LPS and OMP. That is, the primary EcSf2a-2 (cyan) and primary 2457T (magenta) study data from controls who were challenged with wt Shigella are equal to or higher than the corresponding secondary or tertiary exposure data (cyan versus blue, magenta versus red and black) in many cases. We compare modeling results from parameters fit to the primary infection data, which predict theoretical increases in secondary immunity, to those fit to both primary and secondary infection data, which replicate this observed flatness in data. These separate fits enable us to gain insight into the potential effects of establishing LPS and/or OMP immunity at or above measured clinical values and also ensure the robustness of our conclusions across different parameterizations.
It is not known clinically what Ab and/or ASC levels are sufficient for immunity, if any, nor which targets are most immunologically active or involved in protection from disease [3,35]. However, the 27-36% efficacy of the EcSf2a-2 trials indicates that measured levels of ASC, serum IgA, and serum IgG directed against LPS and OMP are not sufficient on average at the group level to prevent infection in the majority of patients. On the other hand, five of six volunteers who received a second wild-type dose of 2457T in the challenge study (that is, 2T: 2457Tx2) were not strongly symptomatic. It is not known if the mechanisms of protection were the measured levels of anti-LPS and anti-OMP serum IgA and/or IgG or other unmeasured immune components perhaps with different bacterial targets. With this mathematical model, we investigate this by using only the measured serum IgA and IgG humoral responses offset to differentiate between trials. Dots represent each individual volunteer's measured values. Means and standard deviations are plotted for each trial group. ET denotes data from the EcSf2a-2 vaccine trials whereas 2T denotes the 2457T rechallenge study. EcSf2a-2 indicates the vaccine was given, 2457T indicates wt challenge, and x1 or x2 indicates once or twice. The five trial groups are the following. ET: 2457Tx1 (cyan): A primary immune response: controls from the EcSf2a-2 trials who were challenged with wild-type (2457T) Shigella without the vaccine. ET: EcSf2a-2x1, 2457Tx1 (blue): An immune response following vaccination: all volunteers who received the EcSf2a-2 vaccine followed by wildtype challenge. 2T: 2457Tx1 (magenta): A primary immune response: unvaccinated volunteers for the 2457T rechallenge study who were challenged once with wild-type 2457T Shigella. 2T: 2457Tx2 (red): A wild-type secondary immune response: volunteers who first served as controls in EcSf2a-2 trials and were challenged with wild-type 2457T (but not vaccinated) in that study. These volunteers were then rechallenged with wild-type 2457T Shigella in the 2457T rechallenge study. 2T: EcSf2a-2x1, 2457Tx2 (black): A secondary immune response after vaccination: volunteers who received both the EcSf2a-2 vaccine and a 2457T challenge in the first studies. They then were rechallenged with wildtype Shigella in the 2457T rechallenge study. They have not been included in the model data fits.
https://doi.org/10.1371/journal.pone.0189571.g001 without any additional immune action and examine in silico whether epithelial infection is sufficiently prevented (indicating hosts with few-to-no symptoms).
To convert clinical data to units more readily usable in the mathematical model, we assume the following. To convert ASC and B M per million peripheral blood mononuclear cells (PBMC) to absolute numbers of cells, we note that there are roughly 10 6 PBMC per mL of whole blood [36] and roughly 5,000 mL of blood in a human body [37]). Thus, we multiply measured ASC or B M cells per million PBMC by a factor of 5,000 to obtain an estimate of cell numbers. We do not have any established standard method to convert our serum IgA and IgG titers to number of antibodies per mL. Thus, we use the following approximation from titer to antibodies per milliliter (Ab/mL) where Ab here is the number of antibodies. IgA and IgG each have a molecular weight of 150 kDa (kg/mol) per Ab [37]. In our data, the measurement of the specific antibodies was quantified and expressed as calculated endpoint titers. However, using a parallel standard curve with known concentration of purified IgG or IgA in the same setup, we approximate that a titer of 1 approximately represents 3 ng/mL of antibody (unpublished data). Taking these factors together with Avogadro's number, we convert our specific antibody titers to the number of specific Ab/mL by a factor of 1: 1.204 × 10 10 Ab/ mL.

Mathematical model
We update our original model in [29] to a system of 12 ordinary differential equations that capture interactions between bacterial and humoral immune dynamics. Biological and mathematical reaction diagrams corresponding to the updated model are given in Fig 2. Model parameter definitions are listed in Table 1. The model focuses on humoral immune agents, specifically Ab, ASC, and B M . We include IgG-and IgA-type immunity for all Ab, ASC, and B M . For simplicity we do not include immunoglobulin M (IgM) production or affinity maturation/somatic hypermutation prior to IgA/IgG secretion, although we consider somatic hypermutation during interpretation of our sensitivity analysis results. We also do not include potential contributions from cytotoxic T-cells. While cell-mediated immunity could help to combat intracellular Shigella within the epithelium, the exact role and importance of effector cell-mediated activity in fighting Shigella infections remains to be well characterized; thus, any mathematical modeling predictions regarding CMI activity against Shigella would be difficult to experimentally corroborate or implement at this time [3]. Our 2013 model [29] used delay differential equations to capture the one-day time delay from naive B cell activation to the production of functional ASC or B M . However, numerical discretization and simulation of these delay differential equations are slow and often stiff due to wide biologically relevant parameter ranges; for numerical efficiency during Monte Carlo simulations for parameter estimation and data fitting, we eliminate the delays and instead use ordinary differential equations. This decision is supported by sensitivity analysis of our original delay model [29], which suggested that these delays have little effect on peak bacterial or antibody magnitudes. By using the best-fit parameter tuples from clinical parameterization of this paper's updated non-delay model to parameterize and simulate a delay differential equation version of our model, we have confirmed that the one-day delayed model produces visually identical simulations to those in this paper (not shown). Despite eliminating these built-in naive cell activation delays, we continue to externally numerically enforce an initial 3.5-day delay after infection begins during which bacteria and any previously established immunity interact without yet eliciting new antigen-induced B cell formation. This time window does not affect pre-existing immunity or innate immune cells such as macrophages, which are modeled indirectly via the engulfed Shigella population.
Our model tracks bacterial populations in the gut lumen, epithelium, and lamina propria after ingestion of contaminated food or water. As an infection begins, Shigella in the lumen, denoted S E , crosses the epithelial barrier at rate μ E via the normal activity of host M cells, which shuttle material from the lumen across the epithelium to the LP. Upon reaching the LP, Shigella, now S I1 , is often engulfed at rate μ NI by nearby innate immune cells such as macrophages that serve as a first line of host defense by engulfing and destroying invaders. However, engulfed Shigella, S N , typically avoids destruction and instead proliferates at rate σ N inside of macrophages. While inside, Shigella is subject to a small rate δ SN of macrophage killing but is safe from antibody targeting. After Shigella breaks out of the macrophage at rate μ N , becoming S I2 , we assume the bacterium is sufficiently distant from M cells and other macrophages to have no likelihood of re-engulfment. We noted in [29] that separate S I1 and S I2 LP populations are mathematically necessary to prevent macrophages from becoming a chronic reservoir for Shigella infections. However, either bacterial population can infect epithelial cells at rate μ I , and once inside, they survive as S C , proliferate at rate σ C , and move freely between epithelial cells without interference from innate or humoral immune defenses. Substantial epithelial cell destruction due to this cellular invasion is primarily responsible for inducing the most severe symptoms of shigellosis [50].
We summarize the bacterial dynamics with the following equations.
The δ terms account for death due to natural causes, spatial washout, or macrophage activity; Shigella removal via antibody is modeled separately.
Interactions between bacteria and the immune system occur when Shigella becomes bound by IgA in the lumen (written A E ) or IgG in the LP (denoted G) at rates α and γ, respectively. Although ASC initially produce IgM before somatically hypermutating and class switching to produce IgA and IgG, for simplicity our model only includes IgA, the most abundant and active antibody isotype at mucosal surfaces, and IgG, the most abundant antibody in serum [37,[51][52][53][54]. We focus on IgA and IgG antibodies that recognize and bind specifically to Shigella outer membrane components, such as LPS or OMP, that are continuously displayed and hence subject to antibody targeting whenever Shigella exists outside of host cells. Equations governing antibody dynamics can be found below. When sufficient numbers (n E or n I , respectively) of IgA or IgG bind to a Shigella bacterium, that bacterium is removed or destroyed. We use mass action interaction terms as a null model to describe bacterial interaction with and removal via antibodies; more complex functional forms are not necessary to match our clinical data. Adding saturating interaction terms would slow immune control of Shigella as antibodies approach a maximum capacity for bacterial neutralization; in comparison, model bacterial dynamics as written could underestimate Shigella numbers and may be a best-case scenario.
The IgG antibody response to Shigella in the LP is largely mediated by IgG-secreting ASC (plasma cells), P G , and upstream IgG-type memory B cells, M G . Similarly, IgA antibodies are created in the LP by IgA-secreting ASC, P A , which can be derived from IgA-type B M cells, M A , or naive B cells (not explicitly modeled); IgA in the LP, A I , then shuttles across the epithelial barrier at rate ω to the lumen where IgA primarily functions, becoming A E . Irrespective of infection state, B M differentiate at baseline rate ϕ 1 to produce ξ 1 ASC, which constantly secrete antibodies at rate β A or β G . However, the presence of bacteria within the lamina propria stimulates naive B cells to create new ASC at rate λ P and new B M cells at rate λ M while also antigenically stimulating additional B M differentiation to create ξ 2 new ASC at rate ϕ 2A or ϕ 2G . The creation of ASC from B M is lessened slightly by Υ and Υ Ã as described below.
These immune dynamics are captured with the following model equations, which together with the bacterial dynamics constitute the full mathematical model.
In the absence of an infection, B M cell generation, death, and differentiation rates must balance to establish an equilibrated immune state with a steady nontrivial (nonzero) memory population. To prevent the mathematical system from becoming neutrally stable at this equilibrium, we assume B M follow logistic growth dynamics with growth rate ρ and carrying capacity κ A or κ G . The carrying capacity is the maximum sustainable population of antigenspecific B M ; that is, it is the size of the Shigella-specific B M population at homeostasis between infections.
In our 2013 paper [29], we found that raising the B M carrying capacity could greatly reduce bacterial load and thus suggested that the B M carrying capacity be investigated further as a potential correlate of Shigella immunity. If B M numbers are a limiting factor in establishing Shigella protective immunity, as indicated by the near-depletion of the B M population in Figs 2 and 3 of that paper [29], sustained B M levels might indeed be crucial. However, the observed depletion could be an artifact of the modeling assumption that B M symmetrically divide and differentiate to form ASC if such differentiation outpaces the assumed logistic growth of B M . Since B M also sometimes act as stem-like cells and divide to replenish their own population [55][56][57], we update the model to include this mixed multiplication strategy. Specifically, all memory differentiation rates are modulated by the percentage Υ or Υ Ã , which are defined as follows. We assume that B M can symmetrically differentiate to produce ξ i ASC (where i = 1 for antigen-independent differentiation and i = 2 for antigen-dependent differentiation) and that this symmetric differentiation occurs Υ% (written as a decimal) of memory differentiations. The remaining 1 − Υ% of differentiations involve a B M asymmetrically dividing to produce one B M replacement and a plasma cell progenitor that immediately differentiates and divides to produce half as many (that is, ξ i /2) ASC. The temporary progenitor cell is not explicitly modeled. Hence, overall a total of Υ% of old B M are lost at the rates given in the model as a total of Υ Ã ξ i ASC per B M are gained, where The degree of depletion observed might depend upon the balance between differentiation and self-replication established by Υ, and thus we examine the robustness of modeling results to a range of values for this parameter.
Because we do not examine cytotoxic T-cell interactions with intracellular Shigella, we also do not model cell-to-cell spread or epithelial cell destruction directly but rather assume high values for S C indicate substantial epithelial cell destruction and severe illness [50]. We also limit our scope to within-host dynamics and do not model any potential bacterial migration from the epithelium back to the lumen (S C ! S E ) before transmitting via a fecal-oral route to another person.

Equilibrium analysis
Bacteria inside the epithelium, S C , are safe from humoral immune targeting and, if present, self-perpetuate and grow without bound in the model, even while the rest of the system reaches equilibrium; a goal of vaccination is to prevent this population from establishing. For equilibrium analysis we leave any growing S C population out through quadrature. The remaining system of equations yields a trivial equilibrium as well as three nontrivial disease-free equilibria: IgA-type only, IgG-type only, and one with both IgA-and IgG-type Ab, ASC, and B M . Nontrivial equilibria represent states in which some level of immunity has been established, although perhaps not at protective levels. We assume both IgA and IgG immune responses are solicited and thus concentrate on the joint IgA/IgG nontrivial equilibrium, for which the equations and parameterized values are given in Table 2.
A positive, stable nontrivial equilibrium is necessary, although perhaps not sufficient, for establishing immunity to Shigella via vaccine or natural infection. Positivity and stability Values of model state variables at the nontrivial disease-free equilibrium for which both IgA and IgG are present are given for the EcSf2a-2 trials or 2457T rechallenge study data. For each, the model is parameterized with data fit to either primary infection data alone or to both primary and secondary infection data. The positive, stable nontrivial equilibrium values are given, representing the long-term presence of some immunity, while in the unstable cases, the trivial equilibrium is instead approached. No positive nontrivial equilibrium was found for parameters that best fit the OMP measurements for the 2457T rechallenge study. The nontrivial equilibrium is evaluated at the corresponding parameter values given in the other tables. The equations for this nontrivial equilibrium are also given, as well as conditions for maintenance of such a nontrivial steady state can be derived by examining the eigenvalues of the Jacobian matrix evaluated at the joint IgA/IgG nontrivial equilibrium. These eigenvalues are given in Table 2. Positivity of the nontrivial equilibrium requires which essentially requires that there be a balance between B M logistic growth and differentiation into ASC; that is, B M must not fully deplete. Stability requires the same condition as positivity plus that This requires that the bacterial growth rate inside of macrophages be less than the rate of bacterial exit so that a chronic or uncontrolled bacterial reservoir is not created inside of macrophages. If stability conditions are satisfied, the nontrivial equilibrium (some immunity) is approached as a stable sink while the trivial equilibrium (no immunity) is a saddle, as are the IgA-only and IgG-only nontrivial equilibria. There is a transcritical bifurcation (not shown) at which the nontrivial IgA/IgG equilibrium becomes negative and switches stability with the trivial equilibrium.

Model implementation
A necessary condition for long-term immunity is a stable nontrivial disease-free equilibrium in the mathematical model. Because the biological threshold for sufficient, i.e. protective, immunity against Shigella is unknown, we explore the protective effects of this nontrivial immune state via our model. First, however, we find parameter value combinations that best fit the model to our clinical data using Latin hypercube sampling and a least squares metric as described below. In the process, we also uncover separate parameter combinations that ensure versus negate stability of the nontrivial immune state; that is, some parameter combinations enable some long-term immunity to establish while other combinations establish no longterm immunity. We numerically investigate the infection dynamics that lead to each. Finally, we use numerical sensitivity analysis to determine which parameter values do and do not control bacterial clearance and the establishment of lasting immunity.
We numerically implement infection dynamics in Matlab using a stiff ODE solver (ode15s). Both primary and secondary infections are simulated so that disease progression from nonimmune versus immune initial states can be compared. We initialize the model using as our initial luminal bacterial level, S E , the Shigella wt or vaccine single-dose levels administered in the clinical trials and reinfections studies: 2 × 10 9 cfu for the EcSf2a-2 vaccine and 1,400 cfu for any wt challenges [30][31][32]34]. Both are above the minimum infectious dose of 100-1,000 cfu [20,31,32,35]. A primary infection is initialized with all other immune and bacterial numbers initially at 0. If a primary infection fails to establish any lasting immunity, as indicated by an unstable nontrivial equilibrium for a given parameterization, we initialize a secondary infection identically to the primary infection. If, on the other hand, the nontrivial IgA/IgG diseasefree equilibrium is stable (with some immunity established), we initialize the secondary infection with immune numbers set to this disease-free equilibrium and then pulse 1,400 cfu of luminal Shigella to start the new infection. Thus, our model initialization mimics the Shigella challenges administered during the clinical studies.
In order to numerically enforce an initial 3.5-day incubation period after infection begins during which bacteria and any previously established immunity interact without yet eliciting new antigen-induced B cell formation, the mathematical model's terms for antigen-dependent B cell proliferation and differentation are tagged and assigned the value zero for the prescribed length of time after initialization of the model simulation, after which these terms switch on to their nonzero values. Such an activation window has no effect on long-term, equilibrium immune levels because bacterial levels are zero at all equilibria and hence there is no chronic antigen induction of immunity. We also investigate how the existence and length of this activation window affect key short-term dynamics and find that, for instance, epithelial Shigella (S C ) levels vary only slightly for between a 0-and 10-day activation delay (not shown).

Parameter estimation
As described earlier, we have serum IgA, serum IgG, IgA-type ASC, and IgG-type ASC time series data from both primary and secondary infections. We focus separately on those that target two bacterial outer membrane components, LPS or OMP, in three different studies, the two EcSf2a-2 vaccine studies and the wild-type Shigella flexneri 2457T challenge study. For the EcSf2a-2 versus 2457T studies separately and for each target separately, our parameter estimations simultaneously fit the model to serum IgA, serum IgG, and both ASC types. It is important to note IgA was not measured in the gut lumen (i.e., stool) of subjects in these studies, and serum IgA is used as an indirect estimate of luminal IgA. We ignore time 0 data points in these fits, as initial measurements cannot always detect low levels of immune components, but we include all other time points. We initially fit only the primary infection with the hope to use the secondary infection data for independent model validation. However, since primary and secondary data from these clinical trials are very similar, as can be seen in Fig 1, the rise in secondary immunity predicted by primary fits does not often match clinical trends in secondary responses. Thus, we additionally fit both primary and secondary data together and compare our results with the primary-only fits to discern general trends between all best fits and to compare immune parameter differences between theorized higher and measured lower secondary immune response values.
In the absence of experimental measurements for many of our rate parameters, we use Monte Carlo parameter estimation runs to extensively explore the multi-dimensional parameter space of unknown parameter values. Of 33 parameters listed in Table 1, we fix seven parameters based on established values in the literature and vary the remaining 26, although some vary together. Tables 3 and 4 give more details. We establish ranges for each parameter, based on literature searches wherever possible, and keep unknown ranges as broad as reason and computation time allow. Using Latin hypercube sampling [58,59], we segment each parameter's range into 10 equal sections and randomly sample a value from each section. We simultaneously sample for all relevant parameters and randomly create a 26-tuple that contains one sampled value for every parameter. To explore a wide range of parameter value combinations, we repeat this 4,000 times. We then run a numerical simulation of the model for both primary and secondary infections to determine model dynamics and data-versusmodel distances for each parameter tuple. We use least squares fitting to determine which parameter tuple best matches our clinical data. In the least squares analysis, we add over the sums of (ln(data) − ln(model)) 2 so that the higher order-of-magnitude antibody values do not mute ASC differences. We save the best-fit tuple and repeat this entire sampling process 25 times to establish a selection of many different parameter combinations that are optimal within their set. This allows us to identify parameter values that are relatively fixed across all best-fit tuples versus parameters that vary more widely without significantly affecting the goodness of fit (S1 Fig in the supporting information). We also look at partial rank correlation coefficients for all of our best-fit cases together to see which, if any, parameters are most crucial in fitting the model to primary and secondary infection data (Fig 3) [60,61]. Finally, a grand Monte Carlo parameter estimation run chooses the single best-fit parameter tuple out of the 1,000,000 total sampled tuples for each study (EcSf2a-2 or 2457T) and each target (LPS or OMP) case.
One pattern we notice is that the best 25 fits in most grand runs generate some best-fit parameter tuples that correspond with a stable positive nontrivial immune equilibrium and others that instead have only the trivial equilibrium where immunity is not lasting. This Parameters values for the model are given along with their ranges and best-fit values using primary infection data only. The model is fit to primary infection data for the EcSf2a-2 trials or the 2457T study using LPS or OMP data separately. For each, the best-fit parameter set is given that produces a stable nontrivial equilibrium (wherein some immunity is established) versus an unstable nontrivial equilibrium (for which no immunity can establish). Parameters with fixed numbers are kept constant (-) while those with values were chosen by fitting from the range provided. Each column should be taken together as the best fit; every individual number is not necessarily the single best-fit value for that parameter. Notes: (1) Lumen Ab: matches the analytical observation that there is a transcritical bifurcation at which the nontrivial equilibrium loses stability. To compare the parameter values and model results associated with a stable versus unstable nontrivial equilibrium, we isolate the best fitting stable nontrivial equilibrium run and the best unstable nontrivial equilibrium run, in which only the trivial Parameters values for the model are given along with their ranges and best-fit values using primary and secondary infection data. The model is fit jointly to both primary and secondary infection data for the EcSf2a-2 trials or the 2457T study using LPS or OMP data separately. For each, the best-fit parameter set is given that produces a stable nontrivial equilibrium (wherein some immunity is established) versus an unstable nontrivial equilibrium (for which no immunity can establish). No parameter set was found that produced a stable nontrivial equilibrium for 2457T OMP data or an unstable nontrivial equilibrium for EcSf2a-2 LPS or OMP data when both primary and secondary infection data are fit. Parameters with fixed numbers are kept constant (-) while those with values were chosen by fitting from the range provided. Each column should be taken together as the best fit; every individual number is not necessarily the single best-fit value for that parameter. Notes: (1) Lumen equilibrium is biologically relevant and stable. This is done for each bacterial target and for each clinical study. The grand best-fit parameter tuples in each case are given in Tables 3 and  4, depending upon if the model is fit to primary infection data only or to both primary and secondary infection data, respectively. Thus, this parameter fitting process gives us a wealth of outputs to parse for identifying the key biological rates and interactions responsible for conferring or failing to confer protective immunity against Shigella.
We use GenSSI to evaluate which parameters are structurally identifiable versus which are difficult to uniquely determine [62]. If all state variables are measurable and if the system launches from a nonimmune state (i.e., the trivial equilibrium), GenSSI finds that δ SE , δ SI , μ E , μ I , μ NI , μ N , ω, β A , β G , α, γ, λ P , and λ M are globally structurally identifiable. No parameters were found to be structurally non-identifiable. However, the identifiability of parameters not listed here could not be resolved due to computational limitations in calculating sufficient Lie derivatives to parse underlying relationships. If we consider only those state variables for which we have data to be measurable, then using 12 Lie derivatives in GenSSI only definitively determines that σ C and δ SC (parameters governing Shigella dynamics inside of epithelial cells) are  Table 1. Extreme PRCC values (> 0.5 or < −0.5) and small p-values indicate parameters that are key for how well the model fits to data. Since all parameters have PRCC values closer to 0 than to ±1, no parameters stand out as most important for these fits. Despite the small PRCC values, some parameters have a p-value less than 0.05: μ I , ξ 2 , ϕ 1 , ϕ 2A , ϕ 2G ; thus, these parameters could have more influence than others in fitting the model to data. Importantly, the parameters most responsible for fitting the model to data could be different than the parameters most responsible for preventing an epithelial infection and resulting symptoms in a host; the latter is explored separately with the sensitivity analysis results.
https://doi.org/10.1371/journal.pone.0189571.g003 structurally non-identifiable. Lack of identifiability is likely for more parameters in this case due to our inability to measure bacterial levels internally in humans.
We conduct numerical sensitivity analysis to discern the effects of individual parameters upon the severity of the Shigella infection, as measured by the number of bacteria within the epithelium, S C , on day 7 after infection. Every parameter is varied one at a time from the model parameterized by the overall best-fit parameter tuple. This is done using as baseline the best-fit tuple for a stable nontrivial equilibrium as well as separately the best-fit tuple for an unstable nontrivial equilibrium, and results are compared when the model fits primary data only versus fits both primary and secondary data. Furthermore, we carry out this sensitivity analysis for each study (EcSf2a-2 or 2457T) and for each bacterial target (LPS or OMP). We compare general trends across all of these scenarios in order to discern which model parameters are key to establishing or preventing severe epithelial infections across all examined biological parameterizations. Parameters are varied at minimum across the ranges given in Table 1 and sometimes across wider ranges. Similar parameters are varied together. The four parameter groups that we identify as sensitive through this analysis under our strictest protection threshold are varied across the following ranges: (0.01, 1)/d for μ E and μ I ; (10 0 , 10 6 ) ASC/ B M on a log scale for ξ 1 and ξ 2 ; (10 2 , 10 10 ) Ab/mL/ASC/d on a log scale for β A and β G ; and (10 2 , 10 8 ) B M on a log scale for κ A and κ G .
The goal of our sensitivity analysis is to identify individual parameters that, when altered, have the capacity to move the resulting dynamics below a protection threshold. Key parameters that do so are potential immune correlates of protection for Shigella and are highlighted for continued investigation in the future. However, it is not known what degree of epithelial invasion by Shigella must be prevented for the host to remain unsymptomatic [3,35]. When estimating this threshold for our model, it should be kept in mind that our model does not incorporate cell-mediated immunity (e.g., T-cells) and thus assumes all bacteria that infect the epithelium escape immune action; in vivo conditions would likely include a naive CMI response that could handle a larger number of bacteria invading the epithelium than the model would suggest. Thus, the protective threshold desired for modeling purposes is the effective number of bacteria in vivo that can not only infect the epithelium by day 7 but also escape all immune action without the host becoming symptomatic. In lieu of a known protective threshold, we examine immune requirements to keep the number of Shigella in the epithelium below either 10 cfu or 100 cfu by day 7 of infection. (Nevertheless, a model threshold of 10 cfu translates to a higher actual number of allowed epithelially invasive bacteria if CMI is present.) Whether immunity that keeps epithelial bacteria below 10 or 100 cfu is protective is not known. Importantly, we will see in our results that the baseline best-fit parameter tuples do not sufficiently prevent the epithelial infection in any case under any parameterization as judged by the 10 cfu threshold. Under the 100 cfu threshold, results will show that three cases would be deemed protective on average with their best-fit parameterizations: (1) fitting to EcSf2a-2 LPS primary immunity data with a stable nontrivial equilibrium, (2) fitting to EcSf2a-2 OMP primary immunity data with a stable nontrivial equilibrium, and (3) fitting to EcSf2a-2 OMP primary and secondary immunity data with a stable nontrivial equilibrium. Since the EcSf2a-2 vaccine demonstrated 27-36% efficacy, it was considered insufficiently protective at the group and population level. This partial efficacy informs the modeling threshold we use as a surrogate for predicting clinical outcome. Since under a 100 cfu threshold the model would have demonstrated complete protection in the three aforementioned cases and since such protection was not demonstrated clinically, we use a 10 cfu threshold as our target level for protection. Nevertheless, we include information using both bacterial thresholds to elucidate how modeling results are impacted by threshold choice.

Primary infection data fits
When we parameterize the model with primary infection data from the vaccine trials (ET: 2457Tx1) or rechallenge study (2T: 2457Tx1), we observe decent fits of the model to primary data, as shown in Panels a and d of Figs 4 and 5 for the EcSf2a-2 vaccine and 2457T rechallenge studies, respectively. This is expected as there are enough unknown parameters that our model has a high likelihood of fitting the data. Some parameter values can range broadly without substantially affecting fits, while other best-fit values cluster together within wider ranges. As seen in S1 Fig, clustered parameter values determine equilibrium stability (σ N and δ SN ), affect bacterial infectivity and Ab migration (μ E , μ I , and ω), and control new B-cell formation (λ P and λ M , which fall broadly in the lower half of their range). We further examine which parameters most affect the model's fit to data by calculating partial rank correlation coefficients (PRCC), which are plotted in Fig 3. When we use PRCC to compare the top 25 best-fit tuples for each case across all cases, no parameter has a PRCC value beyond ±0.5, although some parameters (μ I , ξ 2 , ϕ 1 , ϕ 2A , ϕ 2G ) have a p-value smaller than 0.05. Thus, PRCC analysis indicates that no parameters are the primary contributors to the goodness of fit of the parameterized model to data; a combination of parameters must be responsible for the fits. It is important to note that PRCC analysis describes which parameters are responsible for the model fitting to data, not which parameters are responsible for reducing the epithelial infection or host symptoms; we explore the latter, a key goal of vaccination, with numerical sensitivity analysis.
Our best-fit parameter values when fitting the model to primary infection data are given in Table 3. Given our model's many fitted parameter values, we avoid pitfalls of overfitting by never considering our best-fit parameters to each individually be the best possible parameter value for matching the model to data; rather, we compare commonalities across many best-fit parameter tuples to observe general trends and conclusions that can be made, regardless of whether each parameter value is itself optimized. In this way, we ensure that our conclusions are robust across an array of feasible parameterizations.
We first fit the model to the primary EcSf2a-2 vaccine trial data (ET: 2457Tx1) and separately to the primary 2457T rechallenge data (2T: 2457Tx1). Specifically, we fit the model simultaneously to the serum IgA, serum IgG, A-ASC, and G-ASC data targeting LPS for each study. We separately fit the model to the same data targeting OMP for each study. The fits to primary infection data are plotted in Panels a and d of Figs 4 (EcSf2a-2) and 5 (2457T). We then use the model parameterized against each primary data set to numerically predict bacterial and immune dynamics during a subsequent Shigella infection. The Latin hypercube random samples as a whole were different for all cases, but some of the 1,000,000 tuples generated by Matlab's LHS solver overlapped in parameter values with other cases; surprisingly, the same parameter tuple was optimal for both EcSf2a-2 LPS data and EcSf2a-2 OMP data when fitting the model separately to these primary infection data (but not when fitting to both primary and secondary infection data). This remained true after rerunning these simulations with new random seeds, and thus LPS and OMP EcSf2a-2 primary (but not secondary) data fits have identical characteristics in our results.
The disease-free IgA/IgG nontrivial equilibrium values given in Table 2 are subject to stability conditions, and we analytically and numerically examine whether each best-fit parameterization produces a stable or unstable disease-free state. Model predictions based on fits to primary infection data with predicted secondary infection dynamics are plotted in Panels b, c, e, and f in Figs 4 (EcSf2a-2) and 5 (2457T) and overlaid for comparison (but not fit) with clinical data of secondary immune levels from the two studies. Secondary infection data from the EcSf2a-2 trials show measurements taken after challenging volunteers with wild-type Shigella  EcSf2a-2x1, 2457Tx1); data from the 2457T rechallenge study come from six rechallenge patients who received a sequential double wild-type challenge (2T: 2457Tx2). The EcSf2a-2 studies resulted in 27-36% efficacy at establishing immunity, while five of the six subjects in the 2457T study did not get ill upon rechallenge [32,34].
Results with a stable nontrivial equilibrium when the model is fit only to primary infection data are shown in comparison with measured data values in Panels b and e in Figs 4 (EcSf2a-2) and 5 (2457T). In all cases, the data (taken from Fig 1) have a not-significant difference, and in some cases a decrease, in antibody and immune levels from a primary to a secondary infection. This lack of rise in measured immunity is discordant with the model's predictions in Panel b in Fig 4 and Panels b and e in Fig 5 of heightened secondary immunity arising from a stable nontrivial equilibrium. That is, when a stable nontrivial equilibrium (and hence some immunity) is established by a primary infection, the model overestimates the antibody and ASC levels during a secondary infection in comparison with the data by predicting an increase of roughly 1-3 orders of magnitude in secondary antibody and ASC numbers in all cases except for the EcSf2a-2 OMP response (Panel e in Fig 4). This is consistent with published literature that describes a blunted ASC response to challenge in blood among subjects who are protected from Shigella, presumed to be due to a secondary immune response at the mucosa (as the model predicts) [32]. Interestingly, every grand Monte Carlo fitting round with fits to primary infection data produced a few best-fit parameter tuples that corresponded with an unstable nontrivial equilibrium in addition to those that produced a stable nontrivial equilibrium. Thus, primary infection data can be fit both with parameter combinations that do not lead to an immune state and others that do lead to some immunity. If we look at the parameterizations that produce unstable nontrivial equilibria and initialize our model from the trivial equilibrium (which assumes no lasting immunity from the previous infection), then our model matches the clinical secondary infection data remarkably well for the 2457T data (Panels c and f in Fig 5). This is perhaps unsurprising because the primary and secondary infection data are similar, and the model was fit to the primary data. Nevertheless, Fig 5 reveals that, on average, wt anti-LPS and anti-OMP immune levels during a secondary wt challenge (2T: 2457Tx2) are mathematically nearly identical to initiating an infection from a trivial, non-immune state. That is, the model best matches wt challenge data if no lasting LPS or OMP immunity is established.
The predicted secondary infection levels following a stable disease-free state match data in the case of EcSf2a-2 OMP, which suggests that some anti-OMP immunity (not necessarily protective) was established by the EcSf2a-2 vaccine. Meanwhile, our model predictions fail to match measured secondary infection data for the EcSf2a-2 trials with a stable or unstable equilibrium (that is, with or without some immunity establishing) in every case except when OMP immunity is established (Panel e in Fig 4 versus Panels b, c, and f in Fig 4). The explanation is likely that, although the primary and secondary data are similar, the initial infection dose varies widely in the EcSf2a-2 trials: 2 × 10 9 cfu for a primary vaccine dose versus 1,400 cfu for a secondary challenge dose. Thus, the model did not replicate vaccine trial data when fit solely to primary infection dynamics except when investigating OMP immune responses. We next ask whether we can obtain better information about the EcSf2a-2 vaccine trial by fitting both primary and secondary infection data, and we investigate whether our 2457T and EcSf2a-2 OMP modeling results are robust across other parameterizations. However, our primary data fits inform our primary and secondary data fit results and indicate, for instance, why no best fit parameterizations are found for certain cases (unstable and stable secondary dynamics, respectively) in Figs 6 and 7.

Primary and secondary infection data fits
To see if model results are robust or are artifacts of our parameterization process, we redo the data fits but simultaneously fit both the primary and secondary infection data. We wish to discern if stable nontrivial equilibria, which correspond to establishing some level of immunity (although possibly not protective at the group level), are possible with our model while fitting data other than just EcSf2a-2 OMP, and we seek to corroborate the primary-fit EcSf2a-2 OMP and the primary-fit 2457T modeling results. Thus now in addition to primary infection data (ET: 2457Tx1), we use the data from volunteers who received both the EcSf2a-2 vaccine and a wild-type 2457T Shigella challenge as the secondary infection data (ET: EcSf2a-2x1, 2457Tx1) in the EcSf2a-2 trial fits. For the 2457T secondary infection fits, we match primary infection data (2T: 2457Tx1) as well as data from volunteers who were controls in the EcSf2a-2 trial (and thus were challenged but did not receive the vaccine candidate) and then were rechallenged with wild-type Shigella during the 2457T rechallenge study (2T: 2457Tx2). The results are plotted in Figs 6 (EcSf2a-2) and 7 (2457T) and predicted antibody, ASC, and B M equilibrium numbers are given in Table 2. Best-fit parameter values are listed in Table 4. We again separately examine best-fit parameterizations that produce a stable versus unstable nontrivial equilibrium.
The fits for the OMP measurements from the primary and secondary 2457T study generated only best-fit parameters that corresponded with an unstable nontrivial equilibrium; no stable best fits were found. For every other case, we produce a best-fit parameterization with a stable nontrivial equilibrium that matches both the primary and secondary infection data. In fact, for the EcSf2a-2 trial data (Fig 6), we found no best-fit parameterizations that produced an unstable nontrivial equilibrium and could mimic both primary vaccine data (initialized at 2 × 10 9 cfu) and secondary challenge data (with a 1,400 cfu infection dose). That is, of 25 grand best-fit parameterizations for EcSf2a-2 LPS data and separately for EcSf2a-2 OMP data, all produced only a stable nontrivial immunity. This strongly suggests that the EcSf2a-2 vaccine established some lasting anti-LPS and anti-OMP humoral immunity, although perhaps not at levels that are protective for the population.
Furthermore, we examine the number of Shigella that have invaded the epithelium by day 7 of infection as a proxy for disease severity. This is the day 7 value of the rising purple dashed S C curve in Figs 4-7. Day 7 epithelial Shigella levels are listed in Table 5 for all cases. We note that while all cases have epithelial infections above the 10 cfu threshold, three cases have values below the weaker 100 cfu threshold. These three are (1) fitting to EcSf2a-2 LPS primary infection data with a stable nontrivial equilibrium, (2) fitting to EcSf2a-2 OMP primary infection data with a stable nontrivial equilibrium, and (3) fitting to EcSf2a-2 OMP primary and secondary infection data with a stable nontrivial equilibrium. This too suggests that the EcSf2a-2 vaccine was effective at preventing some disease and establishing some immunity, which is corroborated by the 27-36% vaccine efficacy. However, the level of immunity created by the EcSf2a-2 vaccine was not sufficiently efficacious at a group level. This indicates that even 100 cfu of bacteria (in the absence of a CMI response) invading the epithelium by day 7 are predicted to induce symptoms in the host, so the model's threshold for protection must be lower than 100 cfu of epithelial bacteria. Hence, we use 10 cfu as our target threshold.
Turning to the wild-type 2457T results (Fig 7), the situation is different. When we fit the model to both primary and secondary infection data for a wt challenge, we see that both stable and unstable nontrivial equilibrium cases fit the LPS data well, while only an unstable equilibrium fits OMP data. Furthermore, the nontrivial equilibrium values in Table 2 indicate that the 2457T LPS stable best-fit parameterization (for fitting primary and secondary infection data) establishes negligibly small numbers of Ab, ASC, or B M . Thus, the model only fits the  2457T rechallenge data when little-to-no lasting immunity is established. This is further supported by the day 7 epithelial infection levels (Table 5), which are much larger than 100 cfu for all wt 2457T cases. This would suggest a heavily symptomatic host.
In fact, in all best-fit parameterizations of our model using either EcSf2a-2 or 2457T primary or secondary infection data, bacterial numbers inside of the epithelium rise quickly and uncontrollably. This can be seen by the rising purple dashed S C curves in every panel of Figs 4-7 and by the day 7 levels in Table 5. In three EcSf2a-2 cases previously mentioned, these day 7 epithelial levels are kept below 100 cfu, and yet the vaccine was still not sufficiently efficacious at a group level. In no case is the day 7 epithelial Shigella level kept below the lower threshold of 10 cfu by the humoral immune response. While the true threshold for an asymptomatic host is unknown, these results together indicate that Shigella bacteria are not sufficiently prevented from invading the epithelium by anti-LPS or -OMP immune responses in any investigated scenario or study. Thus, the model predicts epithelial escape of the bacteria (and a symptomatic host) in all best-fit cases regardless of parameterization. That is, in 1,000,000 Monte Carlo runs and 25 best-fit parameter tuples observed for each case, not one predicted that an anti-LPS or anti-OMP vaccine would prevent an epithelial Shigella infection below a 10 cfu threshold, and the three cases that kept the epithelial infection below 100 cfu are known not to have been protective on average at a group level. Thus, our modeling results indicate that humoral immunity against only LPS and OMP, regardless of whether established via vaccine or wild-type infection, will not be protective against Shigella. Below, we amend this slightly while exploring the parameter space further with sensitivity analysis to determine whether altering individual immune response rates can increase the effectiveness of an LPS-or OMP-targeted vaccine.

Sensitivity analysis
To investigate other potential immune targets, we carry out numerical sensitivity analysis on our model and plot key results in Figs 8-11 with corresponding predicted immune levels in S2-S5 Figs. We parameterize our model with the best-fit parameterizations in eight different cases each for EcSf2a-2 and 2457T data: fit to the primary infection only and either producing a stable or unstable nontrivial equilibrium; and fit to both the primary and secondary infection  while the nontrivial equilibrium is either stable or unstable. From each of these best-fit tuples for LPS and separately for OMP, we widely vary one parameter at a time and search for which, if any, can reduce the epithelial Shigella load to below a threshold of 10 cfu on day 7 of a secondary infection. No best-fit parameterizations gave epithelial Shigella values below 10 cfu on day 7, while three EcSf2a-2 best-fit parameterizations (Panels b, h, and i in Figs 8-11) gave values already below a weaker threshold of 100 cfu. Importantly, the true threshold for an asymptomatic host is unknown, so it is possible individuals could be symptomatic with epithelial bacteria below even a 10 cfu threshold; however, establishing a feasible threshold enables us to quantitatively compare cases. We find that only eight (or four groups) of all of our 33 parameters can sufficiently limit the peak epithelial infection to below 10 cfu under one or more parameterizations. Specifically, we see the epithelial infection controlled by drastically reducing the rates that Shigella migrates into the lamina propria (μ I ) or epithelium (μ E ); substantially increasing the rate that B M differentiate into ASC (ξ 1 or ξ 2 ); significantly increasing the rate at which antibodies are produced by activated ASC (β A or β G ); or substantially increasing the Shigella-specific B M carrying capacity (κ A or κ G ), which is the number of Shigella-specific B M at homeostasis between infections. These can be seen in Figs 8-11, which plot peak bacterial levels on day 7 of a secondary infection (following vaccination or infection). The best-fit parameter values in each instance are shown with vertical lines, with actual values given in Tables 3 and 4. We compare the epithelial Shigella curve (solid red) with 10 and 100 cfu thresholds (dashed red lines). The corresponding day 7 peak immune levels of Ab/mL, ASC, and B M are given in S2-S5 Figs. Each group of parameters is varied together; for instance, the listed migration rates (μ terms) are kept identical to one another. All other parameters show little-to-no impact on peak epithelial Shigella numbers and never drop below the 10 cfu epithelial infection threshold despite widely varying their values, sometimes beyond the range of the previous Latin hypercube sampling. Those that drop below the 100 cfu threshold are discussed at the end. Notably, neither Υ (the B M symmetric vs asymmetric division percentage added to this model) nor the length of the incubation period time delay before naive activation affect the degree of epithelial infection (not shown because flat).
The Shigella migration rates to which the model dynamics are sensitive (μ E and μ I ) must be reduced to very small values (Fig 8) to limit the bacterial load within the epithelium; such a reduction would essentially stop the Shigella infection before it starts. However, this would require an immune mechanism different from anti-LPS or -OMP antibody targeting, which is already represented by other terms in the model. Since the causative immune agent is not the currently modeled immune response against LPS or OMP, anti-LPS and -OMP immune levels are nearly flat in almost all cases in S2 Fig despite the potential protection indicated in Fig 8. Notably, the LPS response measured during the wt 2457T rechallenge study is predicted by the model to be insensitive to all parameter variations except a large drop in Shigella migration rates. This indicates that for a wt infection, an anti-LPS antibody response is only effective if it can almost completely block Shigella migration across or into the epithelium (and by a mechanism different than the modeled antibody blocking); otherwise, the model predicts that anti-LPS antibodies are not sufficient for disease prevention by any other mechanism even at much larger Ab/mL, ASC, or B M levels. However, the same is not true for anti-OMP responses or for LPS responses during the EcSf2a-2 vaccine trials.
Epithelial infection levels are sensitive to the ξ 1 , ξ 2 and β A , β G parameter groups in all stable nontrivial anti-OMP cases (Panels h, i, and k in the following: Fig 9 with immune levels in S3  Fig for ξ terms, and Fig 10 and S4 Fig for β terms) and in some but not all anti-LPS responses (Panels b and c in Fig 9 and S3 Fig, and Panel c in Fig 10 and S4 Fig). These parameters are rates at which antibodies and ASC are created, respectively. Oddly, possible protection is observed even in cases without a stable immune state when fitting OMP measurements from the wt 2457T study and varying the β or ξ groups (Panel g in Figs 9 and 10), which suggests large pulses of anti-OMP Ab or ASC production even without vaccination might aid the prevention of symptoms in wt cases. However, this is not observed for LPS or with the EcSf2a-2 vaccine study OMP data. Comparing parameterized values to parameter zones that keep epithelial bacteria below the 10 cfu threshold in all sensitive cases reveals that the production rates of antibodies or ASC must be orders-of-magnitude faster than the best-fit values. Specifically, when ξ 1 and ξ 2 , the rates at which B M produce ASC, are varied, sufficient epithelial protection following the establishment of a stable immune equilibrium occurs with anti-OMP ASC numbers above 10 7 and often greater than 10 9 ASC, with correspondingly high antibody levels at or above 10 15 Ab/mL. However, similar anti-OMP ASC and Ab levels are not protective with an anti-LPS immune response (Panels e and f in Fig 9 and S3 Fig). Varying β A and β G , the rates at which ASC produce Ab, also shows that epithelial protection can occur, but is not guaranteed, when Ab levels rise above approximately 10 15 Ab/mL. Here, a rise in ASC is not required for high Ab levels or protection; in fact, ASC numbers drop in all predicted protective cases as the rate at which ASC produce Ab increases because fewer ASC are required to sustain higher Ab levels.
In some, but not all, best parameter fits shown in Fig 11 and S5 Fig, substantially altering the B M carrying capacities (κ A and κ G ) could reduce or even prevent the epithelial infection; that is, establishing a larger permanent pool of B M could generate protective immunity. This sensitivity is observed when a stable nontrivial equilibrium is established in the EcSf2a-2 vaccine trial as well as for anti-OMP B M in a 2457T wt rechallenge (Panels b, c, h, i, and k in Fig  11). The effect is strongest when the carrying capacity for B M is increased from the best-fit values to establish lasting B M populations of at least 10 6 cells. However, such B M numbers were not reached during infection in the 2457T LPS best-fit case despite raising the B M carrying capacity well above 10 6 cells. This suggests that establishing a larger permanent anti-OMP or anti-LPS B M population above the current levels established with the EcSf2a-2 vaccine candidate or a wt infection can (but does not necessarily) raise overall B M numbers and could correlate with protection in some, but not all, cases.
Comparison of Figs 8-11 in relation to the baseline parameterization of each (vertical lines and Tables 3 and 4) reiterates that all parameterizations that best fit the clinical data have parameters in, and often well into, the zones that produce epithelial infection levels above the 10 cfu threshold, as previously seen in Figs 4-7. If these parameter values can be raised biologically, epithelial invasion by Shigella can be reduced and even potentially entirely prevented. For instance, the epithelial Shigella levels in Panel i of Figs 8-11 drop below 1 cfu if any one of the four sensitive parameter groups is substantially altered. This suggests there might be multiple avenues by which to raise anti-OMP immunity to fully protective levels for the EcSf2a-2 vaccine. In contrast, altering wt 2457T anti-LPS immunity is predicted as never fully protective and wt anti-OMP immunity is fully protective solely when accompanied by a rise in ξ 1 and ξ 2 , the rates at which B M generate ASC. While we change only one parameter group at a time, an unexplored combination of parameter changes might also result in sufficient protection.
If we relax the protection threshold to 100 cfu on day 7, we have noted that three EcSf2a-2 vaccine best-fit parameterizations already are below this. These are the best-fit parameterizations for the EcSf2a-2 study that produced a stable equilibrium when fitting only the primary infection LPS data, only primary infection OMP data, or both primary and secondary infection OMP data. Yet we know the EcSf2a-2 vaccine was not on average protective at a group level, and thus 100 cfu is likely too weak a threshold for protection. Nevertheless, we use this 100 cfu threshold to describe how results change when the threshold considered is raised from 10 cfu. For the 10 cfu threshold, epithelial infection dynamics were sensitive to four parameter groups; for a 100 cfu threshold, the epithelial infection on day 7 can be lowered by changing one of 10 parameter groups. However, this is better parsed by separating the three cases already below threshold (Panels b, h, and i in Figs 8-11) from those which drop newly below threshold given a parameter change. The latter is rare and only occurs for a wt 2457T infection with an unstable equilibrium (no lasting immunity). We see that this case, Panel j in Figs 8-11, has a day 7 epithelial infection that drops below 100 cfu by altering any of the 4 previously sensitive parameter groups or by altering one of 3 new parameter groups as follows (not shown): raising the B M growth rate (ρ), substantially decreasing the percent of B M that act like memory stem cells (Υ), or substantially lowering natural antibody death rates (δ AE , δ AI , and δ G ). It must be noted, however, that altering these parameters does not affect wt 2457T dynamics when a stable nontrivial immune equilibrium is established, and thus this sensitivity is not robust under multiple best-fit parameterizations. It is also interesting to look at which parameters exhibit sensitivity to parameter values for the three cases already below 100 cfu with their best-fit parameterizations (panels bhi). The epithelial Shigella numbers are entirely flat for most parameter values, indicating complete insensitivity of the epithelial infection to those parameters. The parameters to which the epithelial Shigella number show sensitivity in these three cases are the 4 previously sensitive parameter groups plus the following 4 parameter groups (not shown): the antigen-dependent B M differentiation rates, ϕ 2A and ϕ 2G ; the antigen-independent B M differentiation rate, ϕ 1 ; the B M growth rate, ρ; and the neutralization rates of Shigella by Ab, α and γ. The existence of these indicates that if a vaccine can establish sufficient anti-LPS or anti-OMP immunity to drop epithelial Shigella values below 100 cfu, additional immune parameters might gain feasibility as secondary immune targets and as conditional immune correlates.

Discussion
In this paper, we build upon our original mathematical model [29] of the humoral immune response against Shigella by introducing asymmetric division of B M , removing nonessential delays, and comprehensively parameterizing the model with clinical anti-LPS and anti-OMP humoral immune data from three human Shigella studies in the 1990s: two EcSf2a-2 vaccine trials and a subsequent wt 2457T rechallenge study [30][31][32]34]. We use Monte Carlo runs of Latin hypercube sampling with least-squares fitting to parameterize the model both against only primary infection data and against primary and secondary infection data simultaneously for each study and each immune target in order to explore which observed trends are parameterization-dependent versus common across parameterizations. To avoid the pitfalls of overfitting a large number of parameters to data, we do not seek to optimize every parameter value individually but rather consider a large number of best-fit parameter tuples as a whole and identify trends across many best-fit tuples. Our modeling results suggest that, on average, humans would be symptomatic with shigellosis following a purely anti-LPS or -OMP humoral immune response due to an uncontrolled infection of gut epithelial cells that is present across all best-fit model parameterizations. We conduct identifiability analysis to gain insight on which parameter values can be uniquely determined, and we use sensitivity analysis to explore which model parameter values must be altered to prevent epithelial invasion and destruction by Shigella bacteria. We identify four key parameter groups as potential vaccine targets or immune correlates.
Results presented in this paper are only as accurate as the mathematical model, and all models are approximations. Limitations of our approach include utilization of serum IgA as an indirect measure of luminal IgA as well as assessment of immune markers that correlate with protection at the group rather than individual level. Predicted levels of epithelial invasion could be underestimates due to not directly modeling cell-to-cell spread of Shigella within the epithelium, which could exacerbate epithelial bacteria levels even above model predictions. Furthermore, immunological thresholds for protection are not clinically known and thus achieving thresholds set in this paper does not guarantee protection from disease. These limitations are transparently considered in the conclusions made and may be addressed in future studies. Modeling results and predictions should be confirmed with biological studies before being clinically employed. Biological validation of these results is difficult due to limited animal models for Shigella and the necessity for human clinical work. However, a strength of mathematical modeling is the ability to first explore conditions theoretically that are difficult to examine biologically and thereby to inform future biological investigation.
Our best-fit parameterizations that match the EcSf2a-2 vaccine trial data reveal that the vaccine likely established some immunity in some individuals (although not protective at the group level), which makes sense given the 27-36% efficacy of the trials [30][31][32]34]. That is, while the EcSf2a-2 vaccine candidate was not sufficiently efficacious at a group level, the model suggests that measured humoral immune levels during the post-vaccine wt challenge are only possible if some lasting anti-LPS and anti-OMP immunity was established. Nevertheless, the epithelial bacterial numbers (the purple dashed S C curves) in Figs 4 and 6, and Table 5 indicate that hosts vaccinated with an EcSf2a-2 vaccine would still be symptomatic with these established levels of immunity if 100 cfu of bacteria on day 7 is an insufficient threshold for protection (which is suggested by the 27-36% efficacy).
While our modeling suggests the EcSf2a-2 vaccine established some levels of lasting LPS and OMP immunity, the model results also suggest that, to reproduce wt data, little-to-no lasting anti-LPS or -OMP humoral immunity is established by a primary wt 2457T infection. That is, the model matches data best when there is no stable nontrivial equilibrium or when the nontrivial equilibrium has very few Ab, ASC and B M . Since it is known that humans can develop immunity to wt Shigella infections, these modeling results indicate that the primary immune factors responsible for conferring immunity are likely not (or not solely) anti-LPS or anti-OMP humoral immune responses.
Thus, in what is perhaps our strongest modeling result, we observe that in every best-fit model parameterization using any primary and secondary infection data, bacteria invade the epithelium and rise quickly and uncontrollably (Figs 4-7). That is, under the 10 cfu day 7 epithelial infection threshold, the model predicts epithelial escape of the bacteria in every best-fit case regardless of parameterization, bacterial target (LPS or OMP), or clinical study (EcSf2a-2 or 2457T). Since it is the epithelial infection that induces the worst symptoms of shigellosis [50], this uncontrolled epithelial expansion indicates that our model's hypothetical patient (and the average EcSf2a-2 or 2457T study volunteer when protected solely by an anti-LPS or -OMP humoral response) is symptomatic to heavily symptomatic. Thus, our modeling results indicate that a vaccine or natural immunity targeting only LPS and/or OMP will not be protective against Shigella. We amend this slightly after conducting sensitivity analysis to determine whether altering individual immune rates from their parameterized values can increase the effectiveness of an LPS-or OMP-targeted vaccine. Nevertheless, any suggested alterations do not fit observed wild-type or vaccine trial measurements, and it is not yet known if they are biologically feasible.
Putting the ubiquitous rise in epithelial infection levels across all best-fit parameterizations together with the observation that some clinical trial and rechallenge study participants were protected at the measured levels of anti-LPS and -OMP antibodies suggests that healthy individuals were protected at least in part by immune components other than anti-LPS and anti-OMP humoral immune responses. We are unable with these data to determine what those might be, but we use sensitivity analysis to explore possibilities. Our numerical sensitivity analysis reveals that severe epithelial infection is insensitive to all but four parameter groups. That is, limiting the epithelial infection to 10 cfu by day 7 of infection requires radically lowering the rates that Shigella migrates into the lamina propria (μ I ) or epithelium (μ E ); markedly increasing the rate that B M differentiate into ASC (ξ 1 or ξ 2 ); substantially increasing the rate at which antibodies are produced by activated ASC (β A or β G ); or significantly increasing the Shigella-specific B M carrying capacity (κ A or κ G ), that is, the homeostatic Shigella-specific B M population size. Limiting the epithelial infection to below 100 cfu is less difficult, and we also explore the broader set of parameter targets by which to do this; however, those parameter alterations might not be sufficiently protective. As a caveat, the 10 cfu and 100 cfu thresholds for the epithelial infection are the numbers of bacteria in the epithelium that escape all immune action; thus, in future models that include CMI, collaboration between humoral and cellular compartments could provide protection even with more bacteria present.
Decreasing the migration of Shigella from the lumen (that is, lowering μ E ) to prevent Shigella infections from beginning to establish would require, for instance, secretory IgA antibodies that block Shigella entry into the lamina propria via M cells. However, decreased migration from the lumen would require a mode of action by immune components that is different than the modeled mass-action removal of bacteria by LPS-or OMP-directed IgA antibodies, which is not sufficient to prevent infection. On the other hand, lowering the migration rate of Shigella from the lamina propria into the epithelium (that is, μ I ) might be accomplished with antibodies that specifically target the epithelial entry proteins displayed by Shigella for brief periods of time. Such proteins include IpaB/C/D/A, IpgC/D/E, and other known bacterial products [63]. We hope to explore this more in future modeling work that takes into account the short exposure times and brief window for antibody targeting of these proteins. Additionally, measuring secretory IgA in future studies rather than utilizing serum IgA as an indirect estimate will help clarify this issue.
Two parameter groups (ξ 1 , ξ 2 and β A , β G ) to which the epithelial infection levels are sensitive in most stable nontrivial cases control the rates at which antibodies and ASC are created. Orders-of-magnitude faster creation of antibodies or their secreting cells is required for preventing an epithelial infection. This faster production rate maintains total antibody levels at those observed in the data starting from often small pools of B M . It should be noted that we change only one parameter group at a time, and thus smaller increases in combinations of B M , ASC, and antibodies might have similar effects; this possibility was partially explored via the parameter fitting process and did not produce a best-fitting parameter combination within the tuples investigated. The observation from Panel j in Figs 9 and 10 that a substantial decrease in epithelial bacteria might be possible even without a stable immune state by varying the β or ξ groups suggests that mimicking high production rates by prophylactically dosing individuals with large amounts of antibodies targeting OMP, even without creating long-term immunity, could be effective. However, this was observed only when fitting OMP measurements from the 2457T rechallenge study and thus is heavily study-and parameterization-dependent. For vaccine design purposes, more attention should perhaps be paid to the stable nontrivial equilibrium cases that suggest orders-of-magnitude increases in the rates that antibodies (both IgA and IgG) or ASC are produced above all current best-fit levels would be required to keep the epithelial infection below 10 cfu on day 7. It is not known if lowering the epithelial Shigella peak below this dose will be sufficient for protection.
Significantly raising the B M carrying capacities (κ A and κ G ) could limit the epithelial infection in some, but not all, cases. This suggests that establishing a larger permanent pool of B M via vaccination could prevent a severe epithelial infection, if this boosting sufficiently raises the B M carrying capacity in addition to the transitory number of B M . In fact, Panels c and i in Fig  11 predict that the epithelial infection could be almost entirely prevented if the long-term levels of B M directed at OMP could be raised by 2-4 orders of magnitude; this is purely theoretical and remains to be explored experimentally. However, the lack of sensitivity of epithelial bacterial numbers to the LPS B M carrying capacity in the 2457T rechallenge study (Panels d-f in Fig  11) suggests that establishing a larger B M permanent pool will not be sufficient for protection in all cases. Broader protection across all cases is obtained solely by lowering Shigella migration rates into the lamina propria and epithelium.
Sensitivity analysis suggests that creating better antibodies through somatic hypermutation and/or affinity maturation might not be protective, since severe epithelial infection dynamics are insensitive to the parameters that relate to antibody effectiveness (α and γ, not shown). Thus, refining the ability of anti-LPS and anti-OMP antibodies to find, be tuned to, and eliminate Shigella is predicted by modeling to not sufficiently protect against epithelial invasion. Furthermore, comparing requirements to meet the 100 cfu threshold versus the 10 cfu threshold for the day 7 epithelial infection reveals some potential immune correlates that are conditional upon establishing EcSf2a-2 vaccine levels of anti-LPS and anti-OMP humoral immunity. That is, while the unconditionally sensitive immune parameters established by the model are the rate that Shigella migrates into the lamina propria or epithelium, the rate that B M differentiate into ASC, the rate at which antibodies are produced by activated ASC, and the Shigella-specific B M carrying capacity, additional conditional targets for enhancing an already partially protective vaccine include the B M differentiation rates, the B M growth rate, and the neutralization rates of Shigella by Ab. This underscores the idea that multiple immune mechanisms and bacterial targets might work in combination to sufficiently protect against Shigella infections.
In moving forward, modeling conclusions should be tested and confirmed clinically. A planned future modeling direction is to examine Shigella immune protection at the individual rather than group level; that is, we hope to utilize our model to discern immunological differences between those individuals who were protected by the EcSf2a-2 vaccine or wild-type challenge versus those who were not. Additionally, we plan to track epithelial cell damage more overtly in the model in order to directly integrate or predict clinical disease outcomes for each individual. Future modeling should also include new immune targets and correlates that become known in the literature. Measuring and modeling CMI and spatially distinct peripheral versus mucosal immune responses may also assist in a better understanding of the protective immune response against Shigella.
Supporting information S1 Fig. The best-fit parameter tuples resulting from every grand run plotted together. These best-fit parameter values result from fitting the model to every combination of the following four characteristics: EcSf2a-2 or 2457T study, LPS or OMP, fitting primary data (1) or fitting both primary and secondary data (2), and producing a stable or unstable nontrivial equilibrium. The parameters are listed on the horizontal axis in the same order as Table 1. The 3.5-day incubation period delay is included. The parameter tuples are listed in Tables 3 and 4. Black lines show the maximum and minimum of the region explored by Latin Hypercube sampling. The entire space was searched but only the best-fit parameter tuple for each scenario is displayed. The colors are as follows. Magenta: stable, fit (1); Cyan: unstable, fit (1); Red: stable, fit (2); Blue: unstable, fit (2). Stable (some immunity) versus unstable (no immunity) parameter values have been slightly offset horizontally to make them easier to distinguish. Each parameter tuple as a whole is optimized, not each individual parameter. (EPS)  Tables 3 and 4 for the model when fit to (a-c) EcSf2a-2 trials LPS data, (d-f) 2457T study LPS data, (g-i) EcSf2a-2 trials OMP data, (j-k) 2457T study OMP data.
Model predictions for peak immune values by day 7 for Ab/mL, ASC, and B M are shown. These are primary (a,d,g,j) and positive, stable secondary (b,e,h,k) infection dynamics fit to primary infection data only or are positive, stable secondary infection dynamics (c,f,i) resulting from fitting both primary and secondary infection data. Results with an unstable nontrivial equilibrium produced when fitting to both primary and secondary infection data are not shown due to strong similarity to the displayed unstable primary infection best-fit results. No best-fit parameter set for OMP measurements from the 2457T rechallenge study produced a positive stable nontrivial equilibrium when fitting both primary and secondary infection data.
Vertical lines indicate best-fit parameter values. Corresponding peak Shigella numbers are given in Fig 8. (EPS) S3 Fig. Predicted peak immune levels when parameters ξ 1 and ξ 2 are varied together on a log scale from 10 0 to 10 6 ASC/B M . These control the rates at which new ASC are generated by B M in the absence or presence of antigenic stimulation, respectively. All other parameters are fixed to the best-fit parameter values given in Tables 3 and 4 for the model when fit to (a-c) EcSf2a-2 trials LPS data, (d-f) 2457T study LPS data, (g-i) EcSf2a-2 trials OMP data, (j-k) 2457T study OMP data. Model predictions for peak Ab/mL, ASC, and B M values by day 7 are shown. These are primary (a,d,g,j) and positive, stable secondary (b,e,h,k) infection dynamics fit to primary infection data only or are positive, stable secondary infection dynamics (c,f,i) resulting from fitting both primary and secondary infection data. Results with an unstable nontrivial equilibrium produced when fitting to both primary and secondary infection data are not shown due to strong similarity to the displayed unstable primary infection best-fit results. No best-fit parameter set for OMP measurements from the 2457T rechallenge study produced a positive stable nontrivial equilibrium when fitting both primary and secondary infection data. Vertical lines indicate best-fit parameter values. 10 x : the number on the x-axis should be used as the exponent of 10 to obtain the true value. Corresponding peak Shigella numbers are given in Fig 9. (EPS) S4 Fig. Predicted peak immune levels when parameters β A and β G are varied together on a log scale from 10 2 to 10 10 (Ab/mL)/ASC/d. These control the rates at which ASC produce IgA and IgG, respectively. All other parameters are fixed to the best-fit parameter values given in Tables 3 and 4 for the model when fit to (a-c) EcSf2a-2 trials LPS data, (d-f) 2457T study LPS data, (g-i) EcSf2a-2 trials OMP data, (j-k) 2457T study OMP data. Model predictions for peak Ab/mL, ASC, and B M values by day 7 are shown. These are primary (a,d,g,j) and positive, stable secondary (b,e,h,k) infection dynamics fit to primary infection data only or are positive, stable secondary infection dynamics (c,f,i) resulting from fitting both primary and secondary infection data. Results with an unstable nontrivial equilibrium produced when fitting to both primary and secondary infection data are not shown due to strong similarity to the displayed unstable primary infection best-fit results. No best-fit parameter set for OMP measurements from the 2457T rechallenge study produced a positive stable nontrivial equilibrium when fitting both primary and secondary infection data. Vertical lines indicate best-fit parameter values. 10 x : the number on the x-axis should be used as the exponent of 10 to obtain the true value. Corresponding peak Shigella numbers are given in Fig 10. (EPS) S5 Fig. Predicted peak immune levels when parameters κ A and κ G are varied together on a log scale from 10 2 to 10 8 B M . These are the carrying capacities of IgA-type and IgG-type memory B cells, respectively. All other parameters are fixed to the best-fit parameter values given in Tables 3 and 4 for the model when fit to (a-c) EcSf2a-2 trials LPS data, (d-f) 2457T study LPS data, (g-i) EcSf2a-2 trials OMP data, (j-k) 2457T study OMP data. Model predictions for peak Ab/mL, ASC, and B M values by day 7 are shown. These are primary (a,d,g,j) and positive, stable secondary (b,e,h,k) infection dynamics fit to primary infection data only or are positive, stable secondary infection dynamics (c,f,i) resulting from fitting both primary and secondary infection data. Results with an unstable nontrivial equilibrium produced when fitting to both primary and secondary infection data are not shown due to strong similarity to the displayed unstable primary infection best-fit results. No best-fit parameter set for OMP measurements from the 2457T rechallenge study produced a positive stable nontrivial equilibrium when fitting both primary and secondary infection data. Vertical lines indicate best-fit parameter values. 10 x : the number on the x-axis should be used as the exponent of 10 to obtain the true value. Corresponding peak Shigella numbers are given in Fig 11. (EPS)