Modeling Microglia Activation and Inflammation-Based Neuroprotectant Strategies During Ischemic Stroke
MModeling Microglia Activation and Inflammation-BasedNeuroprotectant Strategies During Ischemic Stroke
Sara Amato, Andrea Arnold*
Department of Mathematical Sciences, Worcester Polytechnic Institute, Worcester, MA, USA* corresponding author: [email protected]
Abstract
Neural inflammation immediately follows the onset of ischemic stroke. During this process,microglial cells can be activated into two different phenotypes: the M1 phenotype, which canworsen brain injury by producing pro-inflammatory cytokines; or the M2 phenotype, whichcan aid in long term recovery by producing anti-inflammatory cytokines. In this study, weformulate a nonlinear system of differential equations to model the activation of microglia post-ischemic stroke, which includes bidirectional switching between the microglia phenotypes, aswell as the interactions between these cells and the cytokines that they produce. Further, weexplore neuroprotectant-based modeling strategies to suppress the activation of the detrimentalM1 phenotype, while promoting activation of the beneficial M2 phenotype. Through use ofglobal sensitivity techniques, we analyze the effects of the model parameters on the ratio ofM1 to M2 microglia and the total number of activated microglial cells in the system over time.Results demonstrate the significance of bidirectional microglia phenotype switching on the ratioof M1 to M2 microglia, in both the absence and presence of neuroprotectant terms. Simulationsfurther suggest that early inhibition of M1 activation and support of M2 activation leads to adecreased minimum ratio of M1 to M2 microglia and allows for a larger number of M2 than M1cells for a longer time period.
Keywords: neuroinflammation, ischemic stroke, microglia activation, phenotype switching,neuroprotectants, parameter sensitivity.
Stroke is the second leading cause of death worldwide and the fifth leading cause of death in theUnited States, as well as a leading cause of disability [1–4]. Ischemic strokes account for 87% of allstrokes and are caused by a blockage in a blood vessel due to thrombosis or embolism, resulting inoxygen deprivation of the brain [3]. During an ischemic stroke, cell death and damage occur in theaffected brain area called the infarcted core [5–7]. Immediately following the onset of ischemia, thebody naturally responds with inflammation, which can both worsen brain injury and help in longterm recovery.The goal of this study is to develop a mathematical model of the neuroinflammatory processduring ischemic stroke to analyze both the beneficial and detrimental effects of inflammation. Inparticular, we introduce a new coupled system of nonlinear differential equations to model thedynamic interactions between microglia and cytokines, two of the main components involved inneuroinflammation following stroke onset. Neuroinflammation begins with the activation of mi-croglia, a type of neuroglia residing in the central nervous system [6–11]. This activation peaks1 a r X i v : . [ q - b i o . CB ] O c t round two to three days after stroke and persists for several weeks [12, 13]. Microglia maintainhomeostasis of the brain by continuously monitoring their surrounding environment and respondingto pathological signals released by neighboring brain cells [7, 14, 15]. Based on their type, activatedmicroglia produce either anti-inflammatory cytokines or pro-inflammatory cytokines, thereby caus-ing both beneficial and detrimental effects on the brain post-ischemia.Microglia activation is characterized by two phenotypes: M1 and M2. The M1 phenotype(classical activation) is characterized by the secretion of pro-inflammatory cytokines, which can ex-acerbate the inflammatory response and lead to further brain damage. Pro-inflammatory cytokinesinclude tumor necrosis factor alpha (TNF- α ), interleukin 1 beta (IL-1 β ), nitric oxide, interleukin6 (IL-6), and interleukin 12 (IL-12) [6, 10, 16–21]. Microglia can also be activated into the M2phenotype (alternative activiation) and perform crucial roles in limiting inflammation by releasinganti-inflammatory cytokines, including interleukin 4 (IL-4), interleukin 10 (IL-10), and transforminggrowth factor beta (TGF- β ) [6, 10, 16–19, 21–23].M2 microglia have been shown to dominate at the early stages of inflammation, whereas M1microglia activate more slowly and then become the dominant phenotype for the remainder ofthe neuroinflammatory process [9, 17]. Classical and alternative microglia activation is positivelyinfluenced by the presence of pro-inflammatory cytokines and anti-inflammatory cytokines, respec-tively [10, 19, 21, 23, 24]. Experimental studies have shown that anti-inflammatory cytokines inhibitmicroglia activation to the M1 phenotype and, on the other hand, pro-inflammatory cytokines in-hibit microglia activation to the M2 phenotype [6,10,16,17,21]. Further, experimental studies havealso shown that microglia may switch phenotypes from M1 to M2 and vice versa [10,13,16,19,25–27].The switching from the M2 to M1 phenotype has been cited as an area in need of further re-search [15]. Mathematical models for applications other than ischemic stroke have consideredinteractions between microglia phenotypes but have not included the possibility of switching fromthe M2 to the M1 phenotype [18, 24, 28].In this study, we develop a four-compartment model of microglia and cytokine interaction, whichincludes both the M1 and M2 phenotypes, pro-inflammatory and anti-inflammatory cytokines, andbidirectional phenotype switching between M1 and M2 microglia. Previous mathematical modelsstudying intracellular processes of ischemic stroke inflammation have included terms representinggeneral microglia activation, leukocytes, astrocytes, and neurons [5, 11, 29–32]. However, thesemodels do not consider the two microglia phenotypes or phenotype switching, which we include inthis work to analyze both the beneficial and deleterious effects of microglia activation post ischemicstroke. We also study specifically the effects of M2 to M1 phenotype switching, which may lead toincreased cell damage by bolstered production of pro-inflammatory substances in the brain.Mathematical models of neuroinflammation have included the two microglia phenotypes forapplications other than stroke, including traumatic brain injury (TBI), amyotrophic lateral scle-rosis (ALS), hemorrhagic shock, and Alzheimer’s disease [18, 21, 24]. The models for TBI in [24]and for Alzheimer’s disease in [18] include switching from M1 to M2 but do not include switchingfrom the M2 to the M1 phenotype. The model for ALS presented in [21] includes bidirectionalswitching between microglia phenotypes; however, to the authors’ knowledge, there are no currentmodels of neuroinflammation during ischemic stroke which include bidirectional microglia pheno-type switching. Further studies have explored the interactions between cytokines in general neuralinflammation [8,33] but have not included the interactions of microglia producing these substances.Multiple studies have also explored inflammation with macrophages, which behave in a similarmanner to microglia, in applications such as myocardial infarction [28, 34].Despite the widespread impact of ischemic stroke, there are currently only two clinical treat-ment strategies for clot removal. Tissue plasminogen activator (tPA)-induced thrombolysis is theonly FDA-approved medication to restore blood flow in the brain following ischemia. During this2reatment, tPA is intravenously administered to break up the clot within the blood vessel that iscausing the ischemic stroke. This strategy is limited to a small subset of stroke patients due to itsshort treatment window [13, 35–38]. An alternative to thrombolysis drug treatment is thrombec-tomy, a surgical procedure during which a clot retrieval device is used to mechanically remove theblood clot causing the ischemic stroke. Mathematical models for both thrombolysis drug treat-ment and thrombectomy have been developed, including: compartment models to evaluate theeffects of tPA dose on the effectiveness of treatment [36, 37]; predictive models to identify subsetsof patients who would be eligible for thrombolysis [35]; and a model of clot removal for mechan-ical thrombectomy [39]. Both of these treatment strategies increase the risk for hemorrhage postischemic stroke [40–42].A potential new therapeutic strategy may be to target the distinct microglia phenotypes andpromote M2 activation while simultaneously suppressing M1 activation [7,12,13,25,43]. Recent ex-perimental studies have explored the use of neuroprotective substances such as BHDPC, curcumin,miR-124, salidroside (SLDS), glycine, and celastrol to achieve this aim. An in vitro study showedthat BHDPC, a novel neuroprotectant, was able to promote M2 phenotype polarization [44]. In afollow-up study, BHDPC was shown to reduce the amount of M1 microglial cells and enhance theamount of M2 microglia in middle cerebral artery occlusion-induced ischemic brain in mice aftertreatment with the drug [45]. Curcumin was shown to promote M2 microglia activation and inhibitpro-inflammatory responses both in vitro and in vivo in mice [46]. Injection with the microRNAmiR-124 was shown to decrease of ratio of M1 to M2 microglia in a mouse model [47]. IntravenousSLDS injection decreased M1 microglial cells and increased M2 microglial cells post ischemic strokein a mouse model [48]. The amino acid glycine was shown to promote M2 microglial cells invitro and in vivo in Sprague-Dawley rats [49]. Celastrol was shown to decrease pro-inflammatorycytokines in several studies using rodent models [50, 51].Inspired by neuroprotectant strategies, we further modify the model proposed in this work toinclude time-varying terms aiding in the activation of M2 microglia and inhibiting the activationof M1 microglia. We analyze the effects of these neuroprotectant-based terms on the total numberof activated microglia in the system, with specific interest in the ratio of M1 to M2 microglia, fordifferent simulated treatment onset times. Further, by employing global sensitivity techniques, weanalyze the effects of the model parameters on the ratio of M1 to M2 cells and the total activatedmicroglia in both the absence and presence of the neuroprotectant terms. Results emphasize thesignificance of the microglia phenotype switching rates with respect to model sensitivity whenconsidering the ratio of M1 to M2 microglia, while the number of resting microglia and the microgliaactivation and mortality rates are more significant when considering the total activated microglia.The paper is organized as follows. Section 2 describes the coupled system of nonlinear differentialequations derived to model the interactions between the two microglia phenotypes and the pro- andanti-inflammatory cytokines. Section 3 reviews two global sensitivity analysis techniques, Morriselementary effects and the Sobol method, and provides numerical results when these techniques areapplied to the model derived in Section 2. Section 4 details the neuroprotectant-based terms addedto the model to suppress M1 microglia production and bolster M2 microglia production. Thissection also provides computational simulations of the modified model when the neuroprotectantonset times are varied and sensitivity analysis of the modified model. Section 5 features a discussionof the results and future work, and Section 6 gives a summary and conclusions of this work.3 Model Description
In this section we derive a simplified model of neural inflammation post-ischemic stroke, focusingon the interaction between the M1 and M2 microglia phenotypes and pro- and anti-inflammatorycytokines. We assume a constant source of resting microglia, which activates into the M1 or M2phenotypes immediately following the onset of ischemic stroke. This activation is assumed to occurat a constant rate, influenced by the cytokine concentrations. Bidirectional switching can occurbetween the M1 and M2 microglia phenotypes. Pro-inflammatory ( P ) and anti-inflammatory ( A )cytokines are produced by the M1 and M2 microglia, respectively, further influenced by the currentconcentrations of cytokines. Figure 1 gives a schematic representation of the model.Figure 1: Schematic representation of microglia activation. After the onset of stroke, resting mi-croglia are activated into the M1 or M2 phenotype. Activation to the M1 phenotype is positivelyinfluenced by the presence of pro-inflammatory cytokines ( P ) and negatively influenced by anti-inflammatory cytokines ( A ); vice versa for the M2 activation. M1 microglia release detrimentalpro-inflammatory cytokines, and this production is positively influenced by the concentration ofpro-inflammatory cytokines and negatively influenced by the concentration of anti-inflammatorycytokines. Conversely, M2 microglia produce beneficial anti-inflammatory cytokines, and this pro-duction is positively influenced by the concentration of anti-inflammatory cytokines and negativelyinfluenced by the concentration of pro-inflammatory cytokines. Activated microglia may switchbetween the M1 and M2 phenotypes, with the M1 to M2 switching being positively influenced bythe anti-inflammatory cytokines and the M2 to M1 switching being positively influenced by theconcentration of pro-inflammatory cytokines. 4 .1 Microglia The equations describing the change in M1 and M2 microglial cells are as follows: dM dt = a · k M · H ( P ) · ˆ H ( A ) (cid:124) (cid:123)(cid:122) (cid:125) M1 activation − s M → M · H ( A ) · M (cid:124) (cid:123)(cid:122) (cid:125) M1 to M2 switching + s M → M · H ( P ) · M (cid:124) (cid:123)(cid:122) (cid:125) M2 to M1 switching − µ M · M (cid:124) (cid:123)(cid:122) (cid:125) M1 mortality (1) dM dt = a · k M · H ( A ) · ˆ H ( P ) (cid:124) (cid:123)(cid:122) (cid:125) M2 activation + s M → M · H ( A ) · M (cid:124) (cid:123)(cid:122) (cid:125) M1 to M2 switching − s M → M · H ( P ) · M (cid:124) (cid:123)(cid:122) (cid:125) M2 to M1 switching − µ M · M (cid:124) (cid:123)(cid:122) (cid:125) M2 mortality (2)where a is the number of resting microglia, s M → M and s M → M are constant parameters formicroglia phenotype switching, and k M , k M , µ M , and µ M are constant parameters for theactivation and natural mortality of the microglial cells, respectively. The Hill functions H and ˆ H of the cytokines are of the form H ( x ) = x n x x n x + K n x x (3)and ˆ H ( x ) = K N x x K N x x + x N x x (4)where x is taken to be either P or A (representing the pro-inflammatory or anti-inflammatorycytokines, respectively), n x and N x are constant exponents which control the steepness of thecurves, and K x is the half maximal concentration of the respective cytokine. Note that H ( x ) hasthe form of an increasing sigmoidal curve, whereas ˆ H ( x ) is a decreasing sigmoidal curve. SimilarHill functions have been used in previous modeling of cytokines and cells [8, 18, 21, 24, 28, 34]. Thefollowing subsections detail the model terms for microglia activation and phenotype switching.Table 1 lists the descriptions, units, and nominal values of the corresponding model parameters. Microglia activation
Following the onset of ischemic stroke, resting microglia can become polarized to the M1 andM2 phenotypes [6, 10, 17, 19, 21]. We assume that the resting microglia become activated to eachphenotype at a rate that is influenced by the presence of cytokines. Activation of the microgliato the M1 phenotype is positively influenced by the concentration of pro-inflammatory cytokines[10, 19, 21, 24, 28] and negatively influenced by the concentration of anti-inflammatory cytokines[19, 24]. We use the Hill function H ( P ) to represent the saturating promotion of M1 microglia bythe pro-inflammatory cytokines [14, 52]. Likewise, we use the Hill function ˆ H ( A ) to represent thesaturating inhibition of M1 microglia by the anti-inflammatory cytokines. Similar terms are used torepresent the saturating promotion of M2 microglia by anti-inflammatory cytokines [10,19,21,23,24]and the saturating inhibition of M2 microglia by the pro-inflammatory cytokines [17]. Microglia phenotype switching
There is evidence that microglia may switch phenotypes once activated; however, the switching fromthe M2 to M1 phenotype has been cited as an area for further research [15]. The proposed modelincludes bidirectional switching, so that we may analyze the effects of all possible interactions. Tothis end, we assume that once activated, microglia may switch from the M1 phenotype to the M2phenotype [10, 13, 15, 16, 19–21, 24–27] and from the M2 phenotype to the M1 phenotype [10, 13,16, 19, 25–27]. Since switching from M1 to M2 is positivity influenced by the anti-inflammatorycytokines [16, 21], we multiply this term by the Hill function H ( A ). Likewise, since switching from52 to M1 is positively influenced by the pro-inflammatory cytokines [10, 16, 26], we multiply thecorresponding term by H ( P ). The equations describing the concentration changes of pro-inflammatory ( P ) and anti-inflammatory( A ) cytokines are as follows: dPdt = k P · M · H ( P ) · ˆ H ( A ) (cid:124) (cid:123)(cid:122) (cid:125) pro-inflammatory cytokine production − µ P · P (cid:124) (cid:123)(cid:122) (cid:125) natural decay (5) dAdt = k A · M · H ( A ) · ˆ H ( P ) (cid:124) (cid:123)(cid:122) (cid:125) anti-inflammatory cytokine production − µ A · A (cid:124) (cid:123)(cid:122) (cid:125) natural decay (6)where k P , k A , µ P , and µ A are constant parameters related to the production and decay of pro-inflammatory and anti-inflammatory cytokines. The Hill functions H and ˆ H take the same form asin (3) and (4), respectively. The following subsection details the model terms relating to cytokineproduction, and Table 1 lists the relevant parameter descriptions, units, and nominal values. Cytokine production
We assume that pro-inflammatory cytokines are produced by M1 cells at a rate k P , anti-inflammatorycytokines are produced by M2 cells at a rate k A , and their production is influenced by the presenceof both the pro- and anti-inflammatory cytokines in the system [6, 10, 16–19, 21]. In particular,the concentration of pro-inflammatory cytokines supports additional pro-inflammatory cytokinesproduction and suppresses the production of anti-inflammatory cytokines [6, 10, 16, 17, 21], whichwe model through the use of the Hill functions H ( P ) in (5) and ˆ H ( P ) in (6). Further, the presenceof anti-inflammatory cytokines encourages more anti-inflammatory cytokines to be produced andsuppresses the production of pro-inflammatory cytokines [6, 10, 16, 17, 21], which we model throughthe terms H ( A ) in (6) and ˆ H ( A ) in (5). The Hill functions account for the saturating effects ofthese interactions [14, 52]. In summary, the proposed model describing the interactions between the M1 and M2 microgliaphenotypes and pro-inflammatory and anti-inflammatory cytokines comprises Equations (1), (2),(5), and (6). The complete system of coupled nonlinear differential equations is given by dM dt = a · k M · H ( P ) · ˆ H ( A ) − s M → M · H ( A ) · M s M → M · H ( P ) · M − µ M · M dM dt = a · k M · H ( A ) · ˆ H ( P ) + s M → M · H ( A ) · M − s M → M · H ( P ) · M − µ M · M dPdt = k P · M · H ( P ) · ˆ H ( A ) − µ P · PdAdt = k A · M · H ( A ) · ˆ H ( P ) − µ A · A (7)with H and ˆ H defined as in (3) and (4), respectively.Table 1 lists the index, description, nominal value, and unit for each parameter included in Model(7). Nominal parameter values were chosen to obtain model outputs consistent with experimental6ndex Parameter Description Nominal Value Units1 a Number of resting microglia 1000 cells2 k M Activation rate of microglia to M . k M Activation rate of microglia to M . k P Production rate of P . pg/ml hours · cells k A Production rate of A . pg/ml hours · cells n P Hill coefficient for P in H ( P ) 0 . K P Half-maximal concentration of P pgml n A Hill coefficient for A in H ( A ) 0 . K A Half-maximal concentration of A pgml s M → M Rate of M → M . s M → M Rate of M → M . µ M Mortality rate of M . µ M Mortality rate of M . µ A Natural decay rate of A . µ P Natural decay rate of P . N A Hill coefficient for A in ˆ H ( A ) 0 . N P Hill coefficient for P in ˆ H ( P ) 0 . ® programming language. Specifically, ode15s was utilized to compute the numerical solution of Model (7) using the nominal parametersin Table 1 and the initial conditions M M P (0) = 10 pgml , and A (0) = 10 pgml . Figure 2 shows the resulting model output for the numbers of M1 and M2 microgliaand the concentrations of pro-inflammatory and anti-inflammatory cytokines, as well as the ratioof M1 to M2 cells ( M M
2) and the total number of activated microglia ( M M .
34. The minimum ratiois approximately 0 . . Since the behavior of Model (7) is greatly influenced by the choice of values for the 17 modelparameters, we utilize sensitivity analysis techniques in order to study each parameter’s contributionto model output. In particular, we apply two global sensitivity analysis techniques: the Morrismethod of elementary effects, and the Sobol method. Global sensitivity approaches aim to quantifyhow uncertainty and variability in model output can be attributed to uncertainties in the inputs.We summarize Morris elementary effects and Sobol sensitivity analysis below, describing specificallythe application to this work; for more details on these methods, see [54–58].For each method, consider the nonlinear input-output relation y = f ( q ) , q = [ q , . . . , q ] (8)where y is a scalar response variable and each q i is a model parameter whose index i ( i = 1 , . . . , f ( q ) = ˆ M t ; q ) ./M t ; q ) dt (9) f ( q ) = ˆ ( M t ; q ) + M t ; q )) dt (10)Parameters are admitted to vary over a specified space. In this study, each parameter’s admissiblespace is taken to be the interval of 80 − The Morris method of elementary effects quantifies the effect of varying one parameter at a time ona model output. The interval [0 ,
1] is divided into (cid:96) levels, and parameters are randomly sampledfrom these levels. In this work, we let (cid:96) = 100. Parameters are perturbed one at a time by theincrement ∆ = (cid:96) (cid:96) − . (11)For each parameter q i , we generate two vectors which differ only in the i th entry. Parameters arethen rescaled to their admissible parameter space, and elementary effects are computed by d i ( q ) = f ( q + ∆ e i ) − f ( q )∆ (12)where q is the parameter vector and e i is the i th unit vector. This process is repeated for r = 200samples, and the absolute mean µ ∗ i and variance σ i are computed via the formulas µ ∗ i = 1 r r (cid:88) j =1 | d ji | (13)8 a) M1 microglia. (b) M2 microglia.(c) Pro-inflammatory cytokine concentration. (d) Anti-inflammatory cytokine concentration.(e) Ratio of M1 to M2 cells. (f) Total activated microglia. Figure 2: Numerical solution to the model (7) over the time interval [0 ,
72] hours using the parametervalues in Table 1 and initial conditions M M P (0) = 10 pgml , and A (0) = 10 pgml . The plots in (a)-(d) depict the four model states, while (e) and (f) show the ratioof M1 to M2 microglia and the total number of activated microglia, respectively.9nd σ i = 1 r − r (cid:88) j =1 ( d ji − µ i ) , µ i = 1 r r (cid:88) j =1 d ji (14)The absolute mean µ ∗ i in (13) provides an estimate of the absolute value of the average of elementaryeffects over all samples. The variance σ i in (14) measures how far each elementary effect is from themean. Since large variances indicate dependence on neighboring input values, the variance givesan estimate of the combined effects of the interactions of each parameter with other parameters.In this work, we use the absolute mean µ ∗ i to rank the sensitivity of the parameters. Sobol sensitivity analysis is a variance-based method which quantifies how much of the variabilityin the model output can be attributed to each individual parameter or parameter interactions. Toimplement the Sobol method in this work, parameters are randomly sampled from their respectiveparameter space a total of 50 ,
000 times. Half of these samples form the rows of a matrix A , whichhas dimensions 25 , ×
17 (where 17 is the number of parameters), and the other half form therows of a matrix B , a nonidentical 25 , ×
17 matrix. Seventeen additional 25 , ×
17 matrices,denoted as C , . . . , C , are generated such that each C i corresponds to the parameter q i and hasits i th row taken as the i th row of A and its remaining 16 rows taken from B .A scalar model output is generated for each row of the matrices A , B , and C i for all i by firstrunning a forward simulation of the model with parameter values set to the entries in the respectiverow and then computing the response variable (8). This results in scalar response vectors of size1 × ,
000 for each matrix, denoted by y A , y B , and y C i , where y A = f ( A ) , y B = f ( B ) , y C i = f ( C i ) (15)and the output function f is taken to be either (9) or (10), respectively.The first-order sensitivity indices, S i , and total-order sensitivity indices, S T i , are computedusing the formulas S i = var[ E ( Y | q i )]var( Y ) = ( M · y A · y TC i ) − f ( M · y A · y TA ) − f (16)and S T i = 1 − var[ E ( Y | q ∼ i )]var( Y ) = 1 − ( M · y B · y TC i ) − f ( M · y A · y TA ) − f (17)respectively, where f = 1 M M (cid:88) j =1 y jA · M M (cid:88) j =1 y jB (18)and here M = 25 , S T i to achieve an overall sensitivityranking of the parameters. 10 a) f ( q ) = ´ M t ; q ) ./M t ; q ) dt (b) f ( q ) = ´ ( M t ; q ) + M t ; q )) dt Figure 3: Total-order sensitivity index S T i (blue dot) and absolute mean µ ∗ i (green x) for eachparameter q i in Model (7) with respect to the scalar response variable f ( q ). Parameters are labeledusing the indices i ( i = 1 , . . . ,
17) listed in Table 1.
Figure 3 shows the resulting parameter sensitivity rankings using both the Morris and Sobol meth-ods for the scalar responses given in (9) and (10). Parameters are ranked according to their Morrisabsolute means, µ ∗ i in (13), and total-order Sobol sensitivity indices, S T i in (17). Note that thesensitivity rankings are consistent between methods but depend on the response considered.As shown in Figure 3a, when considering the ratio of M1 to M2 microglia in (9), the microgliaphenotype switching rates s M → M and s M → M are the most sensitive parameters with respect toboth Morris and Sobol sensitivity measures. These are followed by the half-maximal concentrations K A and K P of the anti-inflammatory and pro-inflammatory cytokines, respectively. When insteadconsidering the total activated microglia in (10), Figure 3b shows that most of the parameters havea total-order sensitivity index and an absolute mean very close to 0. The most sensitive parameterfor this response is the number of resting microglia a , followed by the microglia activation rates k M and k M and the microglia mortality rates µ M and µ M . Promoting the activation of M2 microglia while simultaneously suppressing M1 microglia activationhas been cited as a possible neuroprotectant strategy to aid during ischemic stroke [7, 12, 13, 25, 43]and has been the subject of several recent experimental studies [44–51]. To simulate the effectsof such neuroprotectant-based strategies on microglia activation, we modify Model (7) to includetime-varying terms to inhibit the production of M1 microglia and bolster the production of M2microglia cells.In introducing these terms, we focus on analyzing the ratio of M1 to M2 microglia over a 72-hour window post stroke, which has been cited as an important time frame for treatment [10, 12,13, 16, 25, 47, 59, 60]. Once being administered at a specified time during this window, we assumethat the neuroprotectant will have a saturating effect on the activation of M1 and M2 microglialcells for 24 hours, after which the effects will diminish and activation will return to pre-treatmentlevels. We further assume that the neuroprotectant should be administered within 15 hours post11ndex Parameter Description Nominal Value Units18 b Minimum value of N . τ Steepness of N τ Steepness of N L Maximum value of N . t Time at which treatment is applied 0 −
15 hoursTable 2: Indices, descriptions, nominal values, and units of the constant parameters in the time-varying neuroprotectant terms N ( t ) and N ( t ) given in (19) and (20), respectively.stroke in order to extend the early dominance of the M2 phenotype over M1 seen in Model (7) foras long as possible before the number of M1 cells again dominates.Based on these assumptions, we include the following two time-dependent terms: N ( t ) = (cid:40) b + − b e τ t − ( t if t ≤
24 + t b + − b e − τ t − ( t − if t >
24 + t (19)which acts to inhibit M1 activation, and N ( t ) = (cid:40) L e − τ t − ( t if t ≤
24 + t L e τ t − ( t − if t >
24 + t (20)which acts to promote M2 activation. Each term is a continuous piecewise sigmoidal function, where L , b , τ , and τ are constant parameters which control the shape of the respective sigmoid graphs,and t is the onset time of simulated neuroprotectant-based treatment. Note that we account fora time delay of 5 hours in each sigmoid curve reaching its respective point of steepest decline orincline, and we shift the curves when t >
24 + t by 38 hours in order to maintain continuity. Theseterms enter Model (7) as multiplicative factors, with N ( t ) multiplying the M1 activation term in(1) and N ( t ) multiplying the M2 activation term in (2). Corresponding parameter descriptionsand nominal values are given in Table 2.To inhibit the production of M1 microglial cells, N ( t ) has the form of a decreasing sigmoid curveuntil 24 hours after the onset of treatment, where we assume that the effects of the neuroprotectantstart to wear off. After 24 hours post-treatment, N ( t ) becomes an increasing sigmoid curve untilthe M1 activation returns to pre-treatment level. The constant parameter b is the minimum valueof N ( t ) and represents how effective N ( t ) is in suppressing the activation of M1 cells. Note thatif b = 0, N ( t ) would completely turn off the activation of resting microglia to the M1 phenotype.Figure 4a shows N ( t ) for onset times t = 0, 5, 10, and 15 hours.Conversely, to bolster the production of M2 microglia cells, N ( t ) has the form of an increasingsigmoid curve until 24 hours after the onset of treatment, at which point we assume that the effectsstart to wear off. After 24 hours post-treatment, N ( t ) transitions to a decreasing sigmoid curveuntil the M2 activation returns to pre-treatment level. The constant parameter L is the maximumvalue of N ( t ) and represents how effective N ( t ) is in bolstering the activation of M2 cells. Notethat if L = 1, N ( t ) would double the activation of resting microglia to the M2 phenotype. Figure 4bplots N ( t ) when t = 0, 5, 10, and 15 hours. 12 .1 Modified Model Summary and Simulation Results In summary, the modified model for suppressing the activation of the M1 phenotype and bolsteringthe activation of M2 microglia phenotype is given by the following system of equations: dM dt = N ( t ) · a · k M · H ( P ) · ˆ H ( A ) − s M → M · H ( A ) · M s M → M · H ( P ) · M − µ M · M dM dt = N ( t ) · a · k M · H ( A ) · ˆ H ( P ) + s M → M · H ( A ) · M − s M → M · H ( P ) · M − µ M · M dPdt = k P · M · H ( P ) · ˆ H ( A ) − µ P · PdAdt = k A · M · H ( A ) · ˆ H ( P ) − µ A · A (21)Figure 4 plots the neuroprotectant terms N ( t ) and N ( t ) over 72 hours when the simulatedtreatment is administered at t = 0, 5, 10, and 15 hours, along with the corresponding plots ofM1 and M2 cells, the ratio of M1 to M2 microglia, and total activated microglia resulting fromModel (21). Note that in the presence of the neuroprotectant terms, the ratio of M1 to M2microglia remains under one until around 60 hours regardless of treatment onset time. Applyingthe neuroprotectant terms right away (i.e., t = 0) results in the lowest ratio of M1 to M2 microglia.This minimum occurs at 11 . .
5, indicating that, at this time, there are aroundtwice as many M2 microglial cells as there are M1. After this minimum is achieved, the ratioof microglial cells increases and ends at 72 hours with the highest ratio of all the onset timesconsidered. For the other onset times shown, a slightly different behavior is observed: The ratiosincrease prior to treatment onset and then decrease less drastically and level off. Ratios thenincrease and, by the end of the 72 hour period, all end up around one on an upward trend.While not shown here, we note that onset times past the 15 hour mark (i.e., t >
15) yielda similar behavior seen in Figure 4e when t = 5, 10, and 15 hours; however, the ratio of M1 toM2 microglial cells no longer stays under one prior to the simulated treatment onset. Instead, theratios when t >
15 increase prior to the onset time, reach a maximum ratio greater than one, andthen decrease and level off. When the effects of the neuroprotectant begin to wear off, the ratiosincrease and end on an upward trend.
We perform a similar sensitivity analysis on the modified model in (21), utilizing both the Morrisand Sobol methods described in Section 3 to quantify how uncertainty and variability in modeloutput can be attributed to uncertainties in the inputs when adding in the neuroprotectant terms.As before, we consider as scalar response variables the ratio of the M1 to M2 microglia as in (9)and the total activated microglia as in (10).Figure 5 shows the resulting parameter sensitivity rankings from the Morris elementary effectsand Sobol sensitivity analysis when t = 0. Similar results hold when t is varied across theadmissible treatment time. In Figure 5a, when considering the ratio of M1 to M2 microglia in (9),the switching rate of M2 to M1 microglia s M → M stands out as being the most sensitive, followedby the opposite switching rate s M → M and half-maximal concentration K P of pro-inflammatorycytokines. The parameters relating to the neuroprotectant terms (indexed 18-22) are much lesssensitive in comparison.In Figure 5b, when considering the total activated microglia in (10), the number of restingmicroglia a remains the most sensitive parameter. The activation rate k M of the M2 microglia13 a) M1 activation inhibitor. (b) M2 activation promoter.(c) M1 microglia. (d) M2 microglia.(e) Ratio of M1 to M2 cells. (f) Total activated microglia. Figure 4: Numerical solutions to the modified model in (21) over [0 ,
72] hours using the nominalparameters values given in Tables 1 and 2 with t varying from 0 to 15 hours.14 a) f ( q ) = ´ M t ; q ) ./M t ; q ) dt (b) f ( q ) = ´ ( M t ; q ) + M t ; q )) dt Figure 5: Total-order sensitivity index S T i (blue dot) and absolute mean µ ∗ i (green x) for eachparameter q i in Model (21) with respect to the scalar response variable f ( q ). Parameters arelabeled using the indices i ( i = 1 , . . . ,
22) listed in Tables 1 and Table 2 with t = 0.becomes the second most sensitive parameter, followed by the natural mortality rate µ M . Note thatout of the neuroprotectant term parameters, the maximum value L of the N ( t ) term supportingM2 activation is the most sensitive. In this work, we develop a system of four coupled nonlinear differential equations describing thedynamics of activated microglia and cytokines during ischemic stroke. In particular, this modelconsiders the switching of activated microglia between the M1 and M2 phenotypes (in both direc-tions) and lumped compartments representing pro-inflammatory and anti-inflammatory cytokines.Inspired by possible neuroprotectant strategies, additional time-dependent terms are included to aidin the activation of beneficial M2 microglia and inhibit the activation of detrimental M1 microglia.Simulations using nominal parameter values show that the model captures experimentally ob-served behavior of M1 and M2 microglial cells and cytokines post ischemic stroke. Numericalresults further emphasize the importance of bidirectional switching between microglia phenotypes,in particular when considering the ratio of M1 to M2 microglia. Global sensitivity analysis resultsindicate that the parameters relating to phenotype switching in both directions are the two mostinfluential parameters in the absence of terms to suppress M1 microglia and bolster M2 microglia.In the presence of these terms, the switching from M2 to M1 stands out as being the most sensitiveparameter. These results indicate that the rate of switching between phenotypes in both directionsis influential on the overall ratio of M1 to M2 microglia in the system.By including terms to suppress the activation of M1 microglia and bolster the activation of M2microglia, the model results in similar ratios of M1 to M2 microglia as observed in experimentalstudies for neuroprotectants. In particular, in numerical simulations using nominal parametervalues, the modified model with neuroprotectant-inspired terms maintains a ratio of M1 to M2 cellsbelow one for around 62 to 68 hours depending on the onset time. This is a significant extension ofM2 cell dominance over results using the baseline model – in the absence of neuroprotectant-basedterms, the ratio remained under one for only the first 17 hours. Further, early onset time of theneuroprotectant terms leads to a decreased minimum ratio of M1 to M2 microglia, which suggests15ossible early reduction in the detrimental effects of neuroinflammation by maintaining a largeramount of M2 cells for a longer time period.When considering the total amount of activated microglia in the system, global sensitivityresults intuitively show that model parameters relating to the M1 and M2 microglial cells are themost sensitive. These include the number of resting microglia and the activation and mortalityrates of the M1 and M2 cells, while model sensitivity with respect to phenotype switching andcytokine production is negligible. Similar results are seen in both the absence and presence of theneuroprotectant terms, however the parameters relating to M2 become more sensitive than thosefor M1 when neuroprotective terms are included.When instead considering the ratio of M1 to M2 microglia, parameters involving phenotypeswitching and cytokines arise as being the most sensitive. In particular, in the absence of the neu-roprotectant terms, parameters relating to the half-maximal concentrations of anti-inflammatoryand pro-inflammatory cytokines follow the switching parameters as the most sensitive, while thenumber of resting microglia is the least sensitive. Similar sensitivity rankings hold in the presence ofthe neuroprotectant terms, while the parameters relating to the neuroprotectant-inspired functionsare notably some of the least sensitive parameters.Since neuroprotectant strategies aim to decrease the ratio of M1 to M2 microglial cells while notnecessarily altering the total number of activated microglia, it is of interest to further study modelterms and parameters relating to the switching between microglia phenotypes, as well as thoserelating to cytokines. The proposed model can be extended to include separate compartments toaccount for interactions between microglia and specific pro-inflammatory cytokines (e.g., TNF- α )and anti-inflammatory cytokines (e.g., IL-4, IL-10).Future work will be performed to estimate model parameters based on experimental data mea-suring the ratio of M1 to M2 microglia in the absence and presence of neuroprotectant treatment,including the use of time-varying parameter estimation techniques to estimate the N ( t ) and N ( t )terms for specific neuroprotectants without assuming the piecewise sigmoidal forms. Further, com-partments for additional intracellular components (including neurons, astrocytes, and macrophages)will be included to study their interactions with the microglial cells and effects on the ratio of M1to M2 microglia. Additional terms may be included to model existing thrombolysis drug treatmentand explore computational simulation of combination strategies with tPA and novel neuroprotec-tants. Neural inflammation is a natural response following the onset of ischemic stroke, propagated by theactivation of microglia. Resting microglia can become classically activated into the M1 phenotypeand produce detrimental substances or alternatively activated into the M2 phenotype and producebeneficial substances. In this study, we formulate a nonlinear system of differential equations tomodel the interactions between the two microglia phenotypes and the cytokines that each phenotypeproduces during inflammatory response. Additionally, we include terms suppressing the activationof M1 microglia and bolstering the production of M2 microglia to simulate possible neuroprotectantstrategies. Model simulations and global sensitivity analysis results highlight the significance ofbidirectional microglia phenotype switching on the ratio of M1 to M2 microglia, in both the absenceand presence of neuroprotectant-inspired terms. Simulation results further demonstrate that earlyonset of terms to inhibit M1 activation and support M2 activation leads to a decreased minimumratio of M1 to M2 microglia and allows the M2 microglia to dominate the number of M1 microgliafor a longer time window. 16 cknowledgments
This work was partially supported by WPI’s Presidential Fellowship (S. Amato) and by the Na-tional Science Foundation under grant number NSF/DMS-1819203 (A. Arnold).
Declaration of Competing Interest:
None.
References [1] W Johnson, O Onuma, M Owolabi, and S Sachdev. Stroke: A global response is needed.
Bulletin of the World Health Organization , 94(9):634–634A, 2016.[2] PK Singh. World Stroke Day 2019. , October 2019.[3] American Stroke Association. About stroke. ,2020.[4] Centers for Disease Control and Prevention. Stroke facts. , 2020.[5] AJH Newton and WW Lytton. Computer modeling of ischemic stroke.
Drug Discovery Today:Disease Models , 19:77–83, 2016.[6] RA Taylor and LH Sansing. Microglial responses after ischemic stroke and intracerebral hem-orrhage.
Clinical and Developmental Immunology , 2013:746068, 2013.[7] MA Yenari, TM Kauppinen, and RA Swanson. Microglial activation in stroke: Therapeutictargets.
Neurotherapeutics , 7(4):378–391, 2010.[8] WD Anderson, HK Makadia, AD Greenhalgh, JS Schwaber, S David, and R Vadigepalli. Com-putational modeling of cytokine signaling in microglia.
Molecular BioSystems , 11(12):3332–3346, 2015.[9] X Hu, P Li, Y Guo, H Wang, RK Leak, S Chen, Y Gao, and J Chen. Microglia/macrophagepolarization dynamics reveal novel mechanism of injury expansion after focal cerebral ischemia.
Stroke , 43(11):3063–3070, 2012.[10] R Orihuela, CA McPherson, and GJ Harry. Microglial M1/M2 polarization and metabolicstates.
British Journal of Pharmacology , 173(4):649–665, 2016.[11] CD Russo, JB Lagaert, G Chapuisat, and MA Dronne. A mathematical model of inflammationduring ischemic stroke. In
ESAIM: Proceedings , volume 30. EDP Sciences, 2010.[12] Y Lee, SR Lee, SS Choi, HG Yeo, KT Chang, and HJ Lee. Therapeutically targeting neu-roinflammation and microglia after acute ischemic stroke.
BioMed Research International ,2014:297241, 2014.[13] R Guruswamy and A ElAli. Complex roles of microglial cells in ischemic stroke pathobiology:New insights and future directions.
International Journal of Molecular Sciences , 18(3):496,2017. 1714] JH Byrne, R Heidelberger, and MN Waxham.
From Molecules to Networks: An Introductionto Cellular and Molecular Neuroscience . Academic Press, 2014.[15] D Boche, VH Perry, and JAR Nicoll. Activation patterns of microglia and their identificationin the human brain.
Neuropathology and Applied Neurobiology , 39(1):3–18, 2013.[16] X Hu, RK Leak, Y Shi, J Suenaga, Y Gao, P Zheng, and J Chen. Microglial and macrophagepolarization—new prospects for brain repair.
Nature Reviews Neurology , 11(1):56–64, 2015.[17] Y Tang and W Le. Differential roles of M1 and M2 microglia in neurodegenerative diseases.
Molecular Neurobiology , 53(2):1181–1194, 2016.[18] W Hao and A Friedman. Mathematical model on Alzheimer’s disease.
BMC Systems Biology ,10(1):108, 2016.[19] Y Nakagawa and K Chiba. Role of microglial M1/M2 polarization in relapse and remission ofpsychiatric disorders and diseases.
Pharmaceuticals , 7(12):1028–1048, 2014.[20] JD Cherry, JA Olschowka, and MK O’Banion. Neuroinflammation and M2 microglia: Thegood, the bad, and the inflamed.
Journal of Neuroinflammation , 11(1):98, 2014.[21] H Shao, Y He, KCP Li, and X Zhou. A system mathematical model of a cell–cell communi-cation network in amyotrophic lateral sclerosis.
Molecular BioSystems , 9(3):398–406, 2013.[22] A Ledeboer, JJP Brev´e, S Poole, FJH Tilders, and AMV Dam. Interleukin-10, interleukin-4,and transforming growth factor- β differentially regulate lipopolysaccharide-induced productionof pro-inflammatory cytokines and nitric oxide in co-cultures of rat astroglial and microglialcells. Glia , 30(2):134–142, 2000.[23] X Liu, J Liu, S Zhao, H Zhang, E Cai, M Cai, X Ji, RK Leak, Y Gao, J Chen, and X Hu.Interleukin-4 is essential for microglia/macrophage M2 polarization and long-term recoveryafter cerebral ischemia.
Stroke , 47(2):498–504, 2016.[24] LE Vaughan, PR Ranganathan, RG Kumar, AK Wagner, and JE Rubin. A mathematicalmodel of neuroinflammation in severe clinical traumatic brain injury.
Journal of Neuroinflam-mation , 15(1):345, 2018.[25] SC Zhao, LS Ma, ZH Chu, H Xu, WQ Wu, and F Liu. Regulation of microglial activation instroke.
Acta Pharmacologica Sinica , 38(4):445–458, 2017.[26] T Tanaka, K Murakami, Y Bando, and S Yoshida. Interferon regulatory factor 7 participatesin the M1-like microglial polarization switch.
Glia , 63(4):595–610, 2015.[27] C Qin, LQ Zhou, XT Ma, ZW Hu, S Yang, M Chen, DB Bosco, LJ Wu, and DS Tian. Dualfunctions of microglia in ischemic stroke.
Neuroscience Bulletin , 35:921–933, 2019.[28] Y Wang, T Yang, Y Ma, GV Halade, J Zhang, ML Lindsey, and YF Jin. Mathematicalmodeling and stability analysis of macrophage activation in left ventricular remodeling post-myocardial infarction.
BMC Genomics , 13:S21, 2012.[29] MA Dronne, JP Boissel, and E Grenier. A mathematical model of ion movements in greymatter during a stroke.
Journal of Theoretical Biology , 240(4):599–615, 2006.1830] P Orlowski, M Chappell, CS Park, V Grau, and S Payne. Modelling of pH dynamics in braincells after stroke.
Interface Focus , 1(3):408–416, 2011.[31] VL Rayz, L Boussel, MT Lawton, G Acevedo-Bolton, L Ge, WL Young, RT Higashida, andD Saloner. Numerical modeling of the flow in intracranial aneurysms: Prediction of regionsprone to thrombus formation.
Annals of Biomedical Engineering , 36(11):1793–1804, 2008.[32] T Lelekov-Boissard, G Chapuisat, JP Boissel, E Grenier, and MA Dronne. Exploration ofbeneficial and deleterious effects of inflammation in stroke: dynamics of inflammation cells.
Philosophical Transactions of the Royal Society A: Mathematical, Physical and EngineeringSciences , 367(1908):4699–4716, 2009.[33] A Torres, T Bentley, J Bartels, J Sarkar, D Barclay, R Namas, G Constantine, R Zamora,JC Puyana, and Y Vodovotz. Mathematical modeling of posthemorrhage inflammation inmice: Studies using a novel, computer-controlled, closed-loop hemorrhage apparatus.
Shock ,32(2):172–178, 2009.[34] H Malek, MM Ebadzadeh, R Safabakhsh, A Razavi, and J Zaringhalam. Dynamics of theHPA axis and inflammatory cytokines: Insights from mathematical modeling.
Computers inBiology and Medicine , 67:1–12, 2015.[35] DM. Kent, HP Selker, R Ruthazer, E Bluhmki, and W Hacke. The stroke-thrombolyticpredictive instrument: A predictive instrument for intravenous thrombolysis in acute ischemicstroke.
Stroke , 37(12):2957–2962, 2006.[36] A Piebalgs, B Gu, D Roi, K Lobotesis, S Thom, and XY Xu. Computational simulations ofthrombolytic therapy in acute ischaemic stroke.
Science Reports , 8:15810, 2018.[37] B Gu, A Piebalgs, Y Huang, C Longstaff, AD Hughes, R Chen, SA Thom, and XY Xu.Mathematical modelling of intravenous thrombolysis in acute ischaemic stroke: Effects of doseregimens on levels of fibrinolytic proteins and clot lysis time.
Pharmaceutics , 11(3):111, 2019.[38] J Minnerup, BA Sutherland, AM Buchan, and C Kleinschnitz. Neuroprotection forstroke: Current status and future perspectives.
International Journal of Molecular Sciences ,13(9):11753–11772, 2012.[39] G Romero, ML Martinez, J F´elez, G Pearce, and ND Perkinson. Applicability of the GPdevice to the Circle of Willis arteries by using a mathematical model. In , pages 48–53. IEEE, 2011.[40] C Motto, A Ciccone, E Aritzu, E Boccardi, CD Grandi, A Piana, and L Candelise. Hemorrhageafter an acute ischemic stroke.
Stroke , 30(4):761–764, 1999.[41] A Shoamanesh, CS Kwok, PA Lim, and OR Benavente. Postthrombolysis intracranial hem-orrhage risk of cerebral microbleeds in acute stroke patients: A systematic review and meta-analysis.
International Journal of Stroke , 8(5):348–356, 2013.[42] U Neuberger, P Kickingereder, S Sch¨onenberger, S Schieber, PA Ringleb, M Bendszus, J Pfaff,and MA M¨ohlenbruch. Risk factors of intracranial hemorrhage after mechanical thrombectomyof anterior circulation ischemic stroke.
Neuroradiology , 61(4):461–469, 2019.[43] MD Ginsberg. Neuroprotection for ischemic stroke: Past, present and future.
Neuropharma-cology , 55(3):363–389, 2008. 1944] C Li, T Chen, H Zhou, Y Feng, MPM Hoi, D Ma, C Zhao, Y Zheng, and SMY Lee. BHDPCis a novel neuroprotectant that provides anti-neuroinflammatory and neuroprotective effectsby inactivating NF- κ B and activating PKA/CREB.
Frontiers in Pharmacology , 9:614, 2018.[45] C Li, Y Bian, Y Feng, F Tang, L Wang, MPM Hoi, D Ma, C Zhao, and SMY Lee. Neuro-protective effects of BHDPC, a novel neuroprotectant, on experimental stroke by modulatingmicroglia polarization.
ACS Chemical Neuroscience , 10(5):2434–2449, 2019.[46] Z Liu, Y Ran, S Huang, S Wen, W Zhang, X Liu, Z Ji, X Geng, X Ji, H Du, RK Leak, andH Xiaoming. Curcumin protects against ischemic stroke by titrating microglia/macrophagepolarization.
Frontiers in Aging Neuroscience , 9:233, 2017.[47] SH Taj, W Kho, M Aswendt, FM Collmann, C Green, J Adamczak, A Tennstaedt, andM Hoehn. Dynamic modulation of microglia/macrophage polarization by miR-124 after focalcerebral ischemia.
Journal of Neuroimmune Pharmacology , 11(4):733–748, 2016.[48] X Liu, S Wen, F Yan, K Liu, L Liu, L Wang, S Zhao, and X Ji. Salidroside provides neuro-protection by modulating microglial polarization after cerebral ischemia.
Journal of Neuroin-flammation , 15(1):1–11, 2018.[49] R Liu, XY Liao, MX Pan, JC Tang, SF Chen, Y Zhang, PX Lu, LJ Lu, YY Zou, XP Qin,LH Bu, and Q Wan. Glycine exhibits neuroprotective effects in ischemic stroke in rats throughthe inhibition of M1 microglial polarization via the NF- κ B p65/Hif-1 α signaling pathway. TheJournal of Immunology , 202(6):1704–1714, 2019.[50] M Jiang, X Liu, D Zhang, Y Wang, X Hu, F Xu, M Jin, F Cao, and L Xu. Celastrol treatmentprotects against acute ischemic stroke-induced brain injury by promoting an IL-33/ST2 axis-mediated microglia/macrophage M2 polarization.
Journal of Neuroinflammation , 15(1):78,2018.[51] Y Li, D He, X Zhang, Z Liu, X Zhang, L Dong, Y Xing, C Wang, H Qiao, C Zhu, and Y Chen.Protective effect of celastrol in rat cerebral ischemia model: Down-regulating p-JNK, p-c-Junand NF- κ B. Brain Research , 1464:8–13, 2012.[52] G Kleiner, A Marcuzzi, V Zanin, L Monasta, and G Zauli. Cytokine levels in the serum ofhealthy subjects.
Mediators of Inflammation , 2013:434010, 2013.[53] C Ferrarese, P Mascarucci, C Zoia, R Cavarretta, M Frigo, B Begni, F Sarinella, L Frattola,and MGD Simoni. Increased cytokine release from peripheral blood cells after acute stroke.
Journal of Cerebral Blood Flow & Metabolism , 19(9):1004–1009, 1999.[54] A Saltelli, S Tarantola, F Campolongo, and M Ratto.
Sensitivity Analysis in Practice: AGuide to Assessing Scientific Models . Wiley Online Library, 2004.[55] MT Wentworth, RC Smith, and HT Banks. Parameter selection and verification techniquesbased on global sensitivity analysis illustrated for an HIV model.
SIAM/ASA Journal onUncertainty Quantification , 4(1):266–297, 2016.[56] RC Smith.
Uncertainty Quantification: Theory, Implementation, and Applications . SIAM,Philadelphia, 2013.[57] CH Olsen, JT Ottesen, RC Smith, and MS Olufsen. Parameter subset selection techniques forproblems in mathematical biology.
Biological Cybernetics , 113(1-2):121–138, 2019.2058] J Wu, R Dhingra, M Gambhir, and JV Remais. Sensitivity analysis of infectious diseasemodels: Methods, advances and their application.
Journal of The Royal Society Interface ,10(86):20121018, 2013.[59] Q Wang, XN Tang, and MA Yenari. The inflammatory response in stroke.
Journal of Neu-roimmunology , 184(1-2):53–68, 2007.[60] AG Dyker and KR Lees. Duration of neuroprotective treatment for ischemic stroke.