A single-cell mathematical model of SARS-CoV-2 induced pyroptosis and the effects of anti-inflammatory intervention
AA single-cell mathematical model of SARS-CoV-2 induced pyroptosis andthe anti-inflammatory response to the drug tranilast
Sara J Hamis, (cid:89) Fiona R Macfarlane. (cid:89) School of Mathematics and Statistics, University of St Andrews, St Andrews, Scotland, UK. (cid:89)
Both authors contributed equally to this work.
This is version 1 of the pre-print article.
Abstract
Pyroptosis is an inflammatory mode of cell death that contributes to the cytokine storm associatedwith severe cases of coronavirus disease 2019 (COVID-19). Central to pyroptosis induced by severeacute respiratory syndrome coronavirus 2 (SARS-CoV-2) is the formation of the NLRP3 inflammasome.Inflammasome formation, and by extension pyroptosis, may be inhibited by certain anti-inflammatorydrugs. One such drug, tranilast, is currently being evaluated as a COVID-19 treatment target in a clinicaltrial. In this study, we present a single-cell mathematical model that captures the formation of the NLRP3inflammasome, pyroptotic cell death and drug-responses to tranilast. The model is formulated in termsof a system of ordinary differential equations (ODEs) that describe the dynamics of proteins involved inpyroptosis. The model demonstrates that tranilast delays the formation of the NLRP3 inflammasome,and thus may alter the mode of cell death from inflammatory (pyroptosis) to non-inflammatory ( e.g., apoptosis).
Key-words
Pyroptosis, COVID-19, NLRP3 Inflammasome, Cytokine Storm, Tranilast, Mathematical model.
SARS-CoV-2 induced pyroptosis and cytokine storms
COVID-19 is a respiratory illness induced by the coronavirus strain SARS-CoV-2 [1]. Most people infectedby the virus experience mild symptoms, however, in cases of severe disease, life-threatening symptoms canmanifest [2]. These divergent disease trajectories have been, largely, attributed to differences in immuneresponses of infected hosts [3]. When cells register the presence of SARS-CoV-2 virions, a multitude ofhost-protective responses are triggered. Infected cells transmit signals that recruit immune cells, such asmonocytes, macrophages and T-cells, to the site of infection. These immune cells act to eliminate the virusfrom the body, and thus they attack infected cells in which virions may replicate. Current research has shownthat, upon active virion replication and release, SARS-CoV-2 can induce pyroptosis in both epithelial cellsand immune cells [3–9]. Pyroptosis is an inflammatory and rapid mode of cell death that is characterised bythe secretion of pro-inflammatory cytokines, cell swelling and, ultimately, membrane rupture resulting in therelease of cytoplasmic contents to the extracellular environment [10, 11]. Furthermore, cytokines released bypyroptosing cells have been shown to induce pyroptosis in neighbouring bystander cells which can registerthe cytokines’ damage associated molecular patterns (DAMPs) [12]. Fundamentally, both the recruitment ofimmune cells and pyroptotic cell death act to protect the host from the virus. Indeed, in hosts with healthyimmune systems, virus-specific T-cells are recruited to the site of infection [3], and cytokine levels are keptcontrolled by negative feedback regulations [13]. However, if the immune system contesting a viral infec-tion is malfunctioning, a wide array of cytokines may be over-produced due to deregulation of the negativeAugust 11, 2020 1/22 a r X i v : . [ q - b i o . S C ] A ug eedback that controls cytokine levels in healthy immune systems [13]. Such an over-compensating immuneresponse may lead to an uncontrolled cytokine activity, commonly referred to as a cytokine storm [3]. Ele-vated levels of both pro-inflammatory and anti-inflammatory cytokines have been observed in severe casesof COVID-19 [2, 9, 13, 14], commonly manifesting with severe symptoms such as pneumonia and acute res-piratory distress syndrome (ARDS), which may lead to multiple organ failure [13]. Hence, cytokine stormsare associated with poor clinical COVID-19 outcomes [2, 4, 6]. Therefore, suppressing the onset of cytokinestorms and the subsequent inflammation, without completely cancelling-out host-protective effects of theimmune system, is one of the suggested treatment strategies explored to combat COVID-19 [13].One approach to suppress cytokine storms is to inhibit pyroptosis whilst leaving other functionalities ofthe immune system untouched. In fact, there exists a number of anti-inflammatory drugs that can inhibitpyroptosis and, through the inhibition of pyroptosis, alter the cell death mode from inflammatory to non-inflammatory [15]. Non-inflammatory modes of cell death includes apoptosis, which is characterised by cellshrinkage with cell membrane integrity maintained throughout cell death, whereas pyroptosis is characterisedby cell swelling and membrane rupture [16,17]. Inhibiting pyroptosis may provide two key clinical advantages.Firstly, this inhibition can suppress the cytokine storm, and the resulting increased inflammation and tissuedamage that it brings. Secondly, it has been shown that tissue factors released upon pyroptosis may initiateblood coagulation cascades, therefore inhibiting pyroptosis may reduce the risk of blood clotting in COVID-19cases [2]. One drug that inhibits pyroptosis is tranilast, an anti-inflammatory drug that is traditionally usedfor treating inflammatory diseases, such as asthma, in Japan and South Korea [15]. Tranilast is currentlybeing re-purposed for treating COVID-19 related symptoms in a clinical trial [2, 18]. The pathway to pyroptosis
In this subsection we describe key aspects of the intracellular pathway that leads to pyroptosis. For amore comprehensive description of the mechanisms driving pyroptosis, we refer the interested reader todetailed reviews [15, 19–22]. The pathway, as considered in the mathematical model described in Section 2,is illustrated in Figure 1.
Figure 1.
A schematic representation of the pathway to pyroptosis, as considered in our mathematicalmodel. Black arrows represent reactions that involve mass transfer between depicted compounds. Bluearrows represent facilitation of compound formation or reactions. Red arrows represent signals that can beturned on or off in the model. Membrane pores (shown in red) induced by GSDMD-N allows for theoutflux of inflammatory cytokines and the influx of extracellular water, causing the cell to swell and thecell membrane to eventually rupture. Abbreviations are listed in Appendix A.2. Chinese Clinical Trial Registry, registration number: ChiCTR2000030002.
August 11, 2020 2/22entral to pyroptosis is the formation of the inflammasome, a multi-protein complex consisting of threemolecular units (a sensor, an adaptor and an executor) that together enable the release of cytokines and theultimate membrane rupture that characterise pyroptosis. It has been shown that the NLRP3 inflammasomecan be induced by SARS-CoV-2 [6, 7]. This inflammasome is named after its sensor molecule NLRP3, i.e. nucleotide-binding and oligomerisation domain (NBD) leucine-rich repeat (LRR)-containing receptors withan N-terminal pyrin domain (PYD) 3. Thus in this study we focus our attention on the inflammasomecomprising the sensor molecule NLRP3, the adaptor molecule ASC, and the executor molecule caspase-1.Here, ASC is the apoptosis-associated speck-like protein containing a caspase activation and recruitmentdomain (CARD) and a PYD domain, while caspase-1 contains a CARD and a catalytic domain (CD). Asis illustrated in Figure 1, the inflammasome forms via homotypic interactions between NLRP3, ASC andpro-caspase-1. More specifically, NLRP3-ASC and ASC-caspase-1 bindings are facilitated by PYD-PYDand CARD-CARD interactions, respectively. Inflammasome formation can be triggered by the presence ofmicrobial pathogens in a cell, or danger signals from other infected cells. Note that, upon formation, onlyone inflammasome is formed per cell [20].At homeostasis, cells do not contain enough NLRP3 to produce an inflammasome base [21]. Instead, NLRP3levels are elevated in cells when toll like receptors (TLRs) sense DAMPs or pathogen associated molecularpatterns (PAMPs). As SARS-CoV-2 is a positive sense RNA virus, it can indeed be detected by TLRs [2].Upon TLRs detecting DAMPs or PAMPs, cytoplasmic nuclear factor kappa-light-chain-enhancer of activatedB cells (NF- κ B) is translocated to the nucleus, initiating the transcription and, by extension, the synthesisof (inactive) NLRP3. The transcription and synthesis of the pro-inflammatory cytokine pro-interleukin-1 β (pro-IL-1 β ) is also regulated by NF- κ B. The transcription-mediated elevation of (inactive) NLRP3 and pro-IL-1 β is often referred to as the priming step of pyroptosis, which is succeeded by the activation step [20].Inactive NLRP3 can become activated upon stimuli from a wide array intracellular events, including alteredcalcium signalling, potassium efflux and the generation of reactive oxygen species (ROS) [15]. In turn, activeNLRP3 molecules can oligomerise and bind together to form a wheel shaped structure, constituting theinflammasome base. Through homotypic binding, ASC molecules can then bind to the NLRP3 inflamma-some base. Thereafter, pro-caspase-1 can bind to ASC, enabling the proximity-induced dimerisation, andthereby the activation, of caspase-1 [23]. Active caspase-1 then acts to cleave the interleukins pro-IL-1 β andpro-interleukin-18 (pro-IL-18) into their respective mature forms, interleukin-1 β (IL-1 β ) and interleukin-18(IL-18). Caspase-1 also cleaves the protein gasdermin D (GSDMD), releasing the active N-terminal domaingasdermin D (GSDMD-N) which can form pores on the cell membrane. These pores enable the outflux ofthe inflamamtory cytokines IL-1 β and IL-18 (in their mature forms) from the cytoplasm to the externalenvironment [24, 25]. The released cytokines can recruit immune cells to the site of infection and initiatepyroptosis in neighbouring cells. Note that the GDSMD-N derived membrane pores also allow for the influx of extracellular material, e.g. water, to the cell. This influx causes cells to swell until their plasma membraneeventually ruptures, releasing cellular material to the extracellular region. Moreover, data suggests that ASC,pro-caspase-1, GSDMD and pro-IL-18 are present in adequate levels for NLRP3 formation and pyroptosisat homeostasis, and thus these proteins are not upregulated by the TLR-stimulated cytoplasm-to-nucleustranslocation of NF- κ B [19, 26].Once the NLRP3 inflammasome is activated, cell lysis can be delayed but not prevented [27, 28]. How-ever, if the formation of the NLRP3 inflammasome is inhibited, a series of intracellular events resultingin non-inflammatory cell death will instead be initiated. There exists multiple drugs that act to inhibitinflammasome formation in order to prevent pyroptosis, a rigorous list of covalent drugs that target theNLRP3 inflammasome can be found in a review by Bertinara et al. [15]. In this study, we include thepharmacodynamical effects of one of these drugs, namely, tranilast. Tranilast inhibits the formation of theNLRP3 inflammasome by covalently binding to the NBD domain of NLRP3 molecules, thus preventingNLRP3-NLRP3 interactions, and the subsequent formation of the NLRP3 inflammasome [29].August 11, 2020 3/22 igure 2.
NLRP3 inflammasome formation. (a) Through a wide array of stimuli, inactive NLRP3 canbecome activated. These active NLRP3 can then bind together to form the inflammasome base. Theanti-inflammatory drug tranilast inhibits this inflammasome base formation. (b) Once the inflammasomebase is formed, ASC can bind to the NLRP3 molecules in the inflammasome base via homotypicPYD-PYD interactions. Thereafter, via CARD-CARD interactions, inactive caspase-1 molecules can bindto NLRP3-bound ASC molecules, and thereby dimerise and activate. Abbreviations used are listed inAppendix A.2.
Mathematical modelling of pyroptosis
The use of mathematical models to describe apoptotic cell death has been well established and reviewed[30, 31]. However, there are significantly fewer mathematical models that describe pyroptotic cell death.Previous modelling works include more implicit descriptions of pyroptosis [32] and descriptions of specificaspects of inflammasome formation [33–35]. Within this work, we explicitly model the intracellular eventsthat drive SARS-CoV-2 induced pyroptosis, from TLRs detecting DAMPs/PAMPs through to the ultimatemembrane rupture. We study this system in the absence, or presence, of the anti-inflammatory drug tranilast.The goal of our model is to capture key aspects of the pyroptosis pathway that are commonly analysed inexperimental studies. Thus, in this mathematical/computational work, we study the time evolution of nuclearNF- κ B, the NLRP3 inflammasome and its components, the pore-forming protein GSDMD, the inflammatorycytokines IL-1 β and IL-18, and cell swelling.August 11, 2020 4/22 Mathematical model
We formulate a system of ODEs describing the dynamics of the key proteins involved in NLRP3 inflam-masome formation and pyroptosis. Specifically, the model includes the dynamics of NF- κ B (Section 2.1),NLRP3 (Section 2.2), ASC (Section 2.3), Caspase-1 (Section 2.4), GSDMD (Section 2.5), IL-1 β and IL-18(Section 2.6). We also consider changes in cell volume over time (Section 2.7). In Section 2.8, we expandthe model to include pharmacodynamic effects of the anti-inflammatory drug tranilast. Figure 1 providesa pictorial overview of the modelled pathway. Variable names are listed in Table 1 in Appendix A.3 and,throughout this work, the bracket notation [ ] denotes compound concentration. Details concerning mod-elling steps and parameters can be found in Appendices A.3 and A.4, respectively. Numerical simulationresults are provided in Section 3. κ B dynamics
When TLRs register DAMPs or PAMPs, cytoplasmic NF- κ B is translocated to the nucleus in order toinitiate the transcription of inactive NLRP3 and pro-IL-1 β . In the case of SARS-CoV-2, for example, theseDAMPs and PAMPs can be induced by intracellular virion replication and release, or by cytokines fromneighbouring, infected cells. In the model, we consider the cytoplasm-to-nucleus translocation process ofNF- κ B to be initiated by a signal S that is turned on in response to TLR-sensed signals. The concentrationsof cytoplasmic and nuclear NF- κ B are here denoted by [NF- κ B c ] and [NF- κ B n ], respectively, and the rateconstants are k a and k b . Once the inflammasome base has formed we assume that the translocation ofNF- κ B into the nucleus is stopped. Therefore, once F = 0 as defined in Equation (5), the translocationinto the nucleus is prevented, however NF- κ B can still return to the cytoplasm. Thus, the rate of change of[NF- κ B c ] and [NF- κ B n ] are described by the equations[NF- κ B c ] dt = − S F k a [NF- κ B c ] + k b [NF- κ B n ] , (1)[NF- κ B n ] dt = S F k a [NF- κ B c ] − k b [NF- κ B n ] , (2)where we consider the total concentration of NF- κ B to be constant over time [36] so that[NF- κ B c ] + [NF- κ B n ] = constant . (3)The DAMP/PAMP initiated TLR-signal S is modelled as a binary on/off function such that S := (cid:40) on , if a DAMP/PAMP induced signal is present,0 off, otherwise. (4) Upon DAMP/PAMP induced TLR-signalling and nuclear NF- κ B translocation, inactive NLRP3, NLRP3 i , istranscribed and subsequently synthesised. As the transcription and synthesis of inactive NLRP3 is promotedby nuclear NF- κ B, we describe the production rate of NLRP3 using a Hill function. Specifically, a Hillfunction with a constant coefficient rate α , Hill coefficient γ NF and half-max concentration-value NF [37]. Inactive NLRP3 can self-activate in response to a secondary signal, S , which can be turned on bymultiple stimuli such as potassium influx, calcium outflux and abnormal reactive oxygen species (ROS).Once activated, NLRP3 a molecules can oligomerise (bind together) to form a wheel-shaped inflammasomebase, modelled by the concentration of oligomerised NLRP3, i.e., NLRP3 o . In the model, we thus considerNLRP3 protein concentrations in three different forms: [NLRP3 i ], [NLRP3 a ] and [NLRP3 o ]. Inactive andactive NLRP3 decays in the model at a rate δ . The forward reactions from the inactive-to-active, and active-to-oligomerised NLRP3 forms are here assumed to be irreversible, and are respectively described using therates k and k . We additionally make the assumption that the oligomerisation of active NLRP3 will ceaseAugust 11, 2020 5/22nce an adequate amount, n , of NLRP3 has oligomerised and formed a base for the inflammasome, thus weintroduce the step-function F , defined as F = F ([NLRP3 o ]) := (cid:40) o ] < n, o ] = n. (5)We now incorporate all the above mechanisms to describe the rate of change of [NLRP3 i ], [NLRP3 a ] and[NLRP3 o ] over time as, d [NLRP3 i ] dt = α [NF- κ B n ] γ NF NF γ NF + [NF- κ B n ] γ NF − S k [NLRP3 i ] − δ [NLRP3 i ] , (6) d [NLRP3 a ] dt = S k [NLRP3 i ] − k F [NLRP3 a ] − δ [NLRP3 a ] , (7) d [NLRP3 o ] dt = k F [NLRP3 a ] , (8)where the binary signal S is set as S := (cid:40) on , if activation signal is present,0 off, otherwise. (9) Once the inflammasome base is formed, that is when F = 1 as defined in Equation (5), free ASC in the cellcan bind to the NLRP3 inflammasome base. Thus in the model, we consider the concentrations of ASC infree and bound form, denoted [ASC f ] and [ASC b ], respectively. The inflammasome-forming binding is hereconsidered to be irreversible and occurs at a rate k , hence the rate of change of [ASC f ] and [ASC b ] are heredescribed by d [ASC f ] dt = − k (cid:0) − F (cid:1) [NLRP3 o ] [ASC f ] , (10) d [ASC b ] dt = k (cid:0) − F (cid:1) [NLRP3 o ] [ASC f ] . (11)Note that, free ASC levels are adequate for pyroptosis at homeostasis [15, 19], and thus the total amount ofASC is constant, [ASC f ] + [ASC b ] = constant . (12) Pro-caspase-1 can bind to inflammasome-bound ASC upon availability and subsequently dimerise into itsactivated, mature form, caspase-1. In the model, concentrations of pro-caspase-1 and caspase-1 are denotedby [pro-C1] and [C1], respectively. The activation occurs at a rate k and is here assumed to be irreversibleso that the rate of change of [pro-C1] and [C1] can be described by d [pro-C1] dt = − k [ASC b ] [pro-C1] , (13) d [ C dt = k [ASC b ] [pro-C1] . (14)Pro-caspase-1 is present in adequate levels in the cell prior to cellular responses to DAMPs or PAMPs [15,19],therefore the total amount of caspase-1, in both pro- and active form is constant over time so that[pro-C1] + [C1] = constant . (15)August 11, 2020 6/22 .5 Gasdermin dynamics GSDMD-N, which is formed as a result of caspase-1 cleavage of GSDMD, produces the pores on the cellmembrane that are central to the pyroptosis process. In the model, we consider the concentrations of bothGSDMD and GSDMD-N, denoted as [GSDMD] and [GSDMD-N], respectively. We describe the caspase-1facilitated cleavage using a Hill function, where γ C1 is the Hill coefficient, and C1 is the half-max [C1]value. We additionally assume a specific rate for [GSDMD] cleavage, α . Therefore, we describe the rate ofchange of [GSDMD] and [GSDMD-N] in the model, through the equations, d [GSDMD] dt = − α [C1] γ C C1 γ C + [C1] γ C [GSDMD] , (16) d [GSDMD-N] dt = α [C1] γ C C1 γ C + [C1] γ C [GSDMD] . (17)Adequate levels of GSDMD required for pyroptosis are available in the cell at homeostasis [26], therefore weinclude a conservation law, [GSDMD] + [GSDMD-N] = constant . (18) Translocation of NF- κ B, from the cytoplasm to the nucleus, induces the transcription and synthesis of pro-IL-1 β . Pro-IL-18, on the other hand, is not synthesised in response to this NF- κ B translocation. When activecaspase-1 is available, the pro-forms of the interleukins can be cleaved into their activated forms. Subse-quently, once membrane pores have been formed in response to GSDMD-N activity, cytoplasmic interleukinsescape to the extracellular region. We here consider the concentrations of IL-1 β /18 in pro-, cytoplasmicand extracellular form, respectively denoted [pro-IL-1 β /18], [IL-1 β c /18 c ] and [IL-1 β e /18 e ]. In the model,pro-IL-1 β and pro-IL-18 are cleaved by caspase-1 at rates α and α , respectively. Once cleaved, cellularIL-1 β c and IL-18 c escape through the GSDMD-N derived pores at the rates k and k , respectively. IL-1 β : The transcription and synthesis of pro-IL-1 β is promoted by [NF- κ B n ], therefore we describe thisprocess using a Hill function with constant coefficient α , Hill coefficient γ NF and the half-max concentrationof NF- κ B, NF . We additionally consider the decay of IL-1 β at the rate δ . We incorporate the abovemechanisms to describe the rate of change of [pro-IL-1 β ], [IL-1 β c ] and [IL-1 β e ] over time as, d [pro-IL-1 β ] dt = α [NF- κ B n ] γ NF NF γ NF + [NF- κ B n ] γ NF − α [C1] γ C C1 γ C + [C1] γ C [pro-IL-1 β ] − δ [pro-IL-1 β ] , (19) d [IL-1 β c ] dt = α [C1] γ C C1 γ C + [C1] γ C [pro-IL-1 β ] − k G [IL-1 β c ] − δ [IL-1 β c ] , (20) d [IL-1 β e ] dt = k G [IL-1 β c ] . (21)Here, G is a function that allows for [GSDMD-N]-dependent cytokine outflux, and water influx, such that, G = G ([GSDMD-N]) := [GSDMD-N][GSDMD]+[GSDMD-N] . (22) IL-18:
The above mechanisms are incorporated to describe the rate of change of [pro-IL-18], [IL-18 c ] and[IL-18 e ] over time as,August 11, 2020 7/22 [pro-IL-18] dt = − α [C1] γ C C1 γ C + [C1] γ C [pro-IL-18] , (23) d [IL-18 c ] dt = α [C1] γ C C1 γ C + [C1] γ C [pro-IL-18] − k G [IL-18 c ] , (24) d [IL-18 e ] dt = k G [IL-18 c ] . (25)The levels of pro-IL-18 in the cell are kept at a homeostatic level [15,19], therefore we include a conservationlaw, [pro-IL-18] + [IL-18 c ] + [IL-18 e ] = constant . (26) The GSDMD-N derived membrane pores that allow for the outflux of mature interleukins, also allow for theinflux of extracellular material ( e.g. water). This influx causes the cell to swell until it eventually ruptures.Single-cell analysis has revealed that before the ultimate membrane rupture occurs, the cell volume increasesgradually [38]. Thus, in the model, we consider the volume of the cell, V , to increase at a rate k once poresare formed on the cell membrane. We describe the rate of change of the cell volume by the equation, dVdt = k G V. (27)Once the cell volume reaches a critical volume V c then the cell ruptures and all cell processes cease. The anti-inflammatory drug tranilast (TR) binds to the NBD domain of active NLRP3 and thus inhibits theinflammasome formation by preventing NLRP3-NLRP3 interactions [29]. Thus in the model, we considerTR binding to NLRP3 a to form the complex TR · NLRP3 a . The forward and reverse rate constants for thisreaction are denoted by k + T R and k − T R respectively. If we include these drug mechanisms in our model,Equation (7) describing the dynamics of NLRP3 a must be modified to include two additional terms, d [NLRP3 a ] dt = . . . − k +TR [TR] [NLRP3 a ] + k − TR [TR · NLRP3 a ] . (28)The rate of change of [TR] and [TR · NLRP3 a ] can now be included in the model as, d [TR] dt = − k +TR [TR] [NLRP3 a ] + k − TR [TR · NLRP3 a ] , (29) d [TR · NLRP3 a ] dt = k +TR [TR] [NLRP3 a ] − k − TR [TR · NLRP3 a ] , (30)where the following conservation law holds,[TR]+[TR · NLRP3 a ] = constant. (31) The model is implemented in
MATLAB [39], where the system of ODEs (1)-(31), given in full in Ap-pendix A.1, are solved numerically using the function ode15s . A more detailed description of the numericalset-up can be found in Appendix A.5, which also includes instructions on how to access and run the code.Simulation results are provided and discussed in the next section.August 11, 2020 8/22
Results and discussion
We first consider the model without any drug, that is, we simulate the system given by Equations (1)-(27).Simulation results are provided in Figure 3, in which compound concentrations (in arbitrary units) and cellvolume are plotted over time. These results show that cytoplasmic NF- κ B is quickly translocated to thenucleus upon DAMP/PAMP-induced TLR signalling as [NF- κ B c ] quickly decreases in favour of [NF- κ B n ]at the start of the simulation. Later on, once the inflammasome base has formed, i.e., when [NLRP3 o ] hasreached the value n , cytoplasmic NF- κ B translocates to the nucleus at a decreased rate, thus halting thetranscription, and by extension the synthetisation, of NLRP3 i and pro-IL-1 β . Synthesised NLRP3 i activatesto become NLRP3 a , which in turn oligomerises and binds together to form the inflammasome base, hereexpressed in terms of the concentration of NLRP3 o . These dynamics are captured in the subplot displayingthe time evolution of [NLRP3 i ], [NLRP3 a ] and [NLRP3 o ]. Note that once the inflammasome base is formed,the NLRP3 a -to-NLRP3 o oligomerisation stops and the binding of ASC to the inflammasome is initiated.This can be observed in the results, as [ASC b ] levels increase when [NLRP3 o ] reaches the value n . Further,as [ASC b ] levels increase, the dimerisation and activation of caspase-1 is facilitated so that [pro-C1] levelsdecrease and [C1] levels increase. In turn, when [C1] levels increase, the cleavage of GSDMD, pro-IL-1 β andpro-IL-18 occurs. This is apparent in the subplots, as [GSDMD-N], [IL-1 β c ] and [IL-18 c ] levels start increas-ing once [C1] levels surpass zero. The outflux of mature interleukins and the influx of extracellular water tothe cell is mediated by GSDMD-N derived pores, here implicitly modelled via the concentration [GSDMD-N].Thus, once GSDMD-N has been cleaved, so that [GSDMD-N] >
0, the cytoplasmic interleukin concentrations[IL-1 β c ] and [IL-18 c ] respectively decrease in favour of the extracellular concentrations [IL-1 β e ] and [IL-18 e ].[GSDMD-N] levels over zero also enable the cellular volume increase that eventually results in membranerupture once the volume reaches the maximal capacity V c . Note that when membrane rupture occurs, allintracellular mechanisms are ceased and thus the time progression in the subplots stops. Figure 3.
Numerical simulations results of the ODE model described by Equations (1)-(27) - that is, themodel of pyroptosis in the absence of drugs. We display the concentration of each model component inarbitrary units (a.u) over time, as well as cell volume dynamics. Note that the time progression stops oncethe critical volume V c is reached and the ultimate cell membrane rupture occurs.August 11, 2020 9/22e next consider numerical simulations of the model including the effects of the drug tranilast, that isthe system given by Equations (1)-(27), with the modifications and additions described in Equations (28)-(30). The results of the numerical simulations are displayed in Figure 4, where we compare different initialconcentrations (dosages) of the drug. We display the time evolution of each model component in a separatesubplot, where the graph colour corresponds to drug dosage. For all components, increasing levels of thedrug push pyroptotic events further in time. Note that the results for the maximal tested drug dosageresults in many of the processes driving pyroptosis not occurring in the time frame examined. The drugtranilast specifically inhibits the NLRP3 a -to-NLRP3 o oligomerisation process. Therefore, when tranilast isadded to the system, we observe that as the drug concentration [TR] increases, so does the amount of timeit takes for [NLRP3 o ] to reach the threshold value n . This means that the inflammasome base formation,and by extension the downstream processes resulting in inflammatory cytokine secretion and the ultimatemembrane rupture, are delayed. Recall that if pyroptosis is delayed in a dying cell, the cell death mode canbe altered from inflammatory to non-inflammatory. Figure 4.
Numerical simulations results of the ODE model including the drug tranilast. We display theconcentration of each model component in arbitrary units (a.u) over time, as well as cell volume dynamicsfor several dosages of the drug. The dynamics of each model component is plotted in a separate subplot,where the line colour corresponds to applied drug dosage, as described in the legend. Note that the timeprogression stops once the critical volume V c is reached and the ultimate cell membrane rupture occurs. Pyroptosis has been identified as a key mechanism involved in the cytokine storm and the elevated inflam-mation associated with severe cases of COVID-19 [2]. Consequently, pyroptosis has been suggested as aCOVID-19 treatment target. In order to investigate the pathway of pyroptosis in further detail, we have inthis study formulated a single-cell mathematical model that captures the key proteins and processes involvedin this pro-inflammatory mode of cell death. Specifically, we model the process from DAMP/PAMP-inducedTLR signalling to cell death and membrane rupture. The model is formulated in terms of a system of ODEsAugust 11, 2020 10/22escribing the dynamics of protein concentrations and cell volume. We further expand the model to includepharmacodynamic effects of the anti-inflammatory drug tranilast, and our numerical simulations demonstratethat increasing drug levels result in delayed NLRP3 inflammasome formation, and by extension, delayed cellrupture. Recently it has been found that if GSDMD-N levels in the cell are low, inflammasome activationcan lead to apoptosis instead of pyroptosis as GSDMD-N inhibits caspase-3 action [11, 19, 22, 40, 41]. There-fore, by delaying the cleavage of GSDMD, the cell death mechanism may be switched from inflammatoryto non-inflammatory. As we have shown in our simulations inclusion of the drug tranilast can delay theformation of the inflammasome, and therefore may provide a method of switching cell death mode.In this paper we have provided mathematical model for intracellular pyroptosis dynamics in a non-specifiedcell subjected to a non-specified pyroptosis-inducing stimuli. Our model parameter values were motivatedby data from a collection of previous studies with different experimental settings. To quantitatively capturepyroptosis in a specific setting, we would require access to comprehensive data to update parameter choicesaccordingly. Furthermore, comprehensive data would also allow us to evaluate absolute concentrations foreach of the proteins included within the model. With experimentally validated parameter values, a mean-ingful sensitivity analysis could be performed. However, as described in the Acknowledgements, this work ispart of the
Rapid Assistance in Modelling the Pandemic (RAMP) initiative and thus we choose to make ourcurrent work available now.In regards to further applications of this work, the single cell model developed in this study could be extendedto investigate other anti-inflammatory drugs that target the NLRP3 inflammasome, a comprehensive reviewof such drugs is provided by Bertinara et al. [15]. This single cell model could also be incorporated intomulti-cell models in order to study the effect of pyroptosis on a cell population or tissue scale. This wouldbe especially interesting in a multi-scale model incorporating both epithelial cells and immune cells, as therelease of cytoplasmic cytokines such as IL-1 β and IL-18 into the external environment has been linkedto contribute to the cytokine storm. It has been shown that IL-1 β can induce pyroptosis in neighbouringcells [12] and that IL-18 can induce interferon (IFN) γ production in neighbouring cells, which plays a role inthe recruitment of immune cells [42]. Therefore, it would be interesting to study these bystander effects in amulti-cell setting. Accordingly, we are collaborating with an interdisciplinary coalition of mathematicians andbiologists to model SARS-CoV-2 induced within-host dynamics [43]. This multi-scale model is implementedin the PhysiCell [44] tissue simulator, allowing for in silico investigations of COVID-19 related dynamics ofepithelial cells, immune cells, virions and cytokines. Acknowledgements
The authors would like to thank the SARS-CoV-2 Tissue Simulation Coalition for ongoing collaboration andfeedback on model development. Further, this work is part of the RAMP (Rapid Assistance in Modelling thePandemic) initiative, coordinated by the Royal Society, UK. The authors would like to thank Prof. Mark A.J.Chaplain (University of St Andrews) for coordinating our RAMP Task Team modelling within-host dynamics.
SARS-CoV-2 Tissue Simulation Coalition: http://physicell.org/covid19/ .RAMP: https://royalsociety.org/topics-policy/health-and-wellbeing/ramp/ . Funding
SJH was supported by the Medical Research Council [grant code MR/R017506/1].August 11, 2020 11/22 eferences [1] P. Zhou, X. Yang, X. Wang, B. Hu, L. Zhang, W. Zhang, H. Si, Y. Zhu, B. Li, C. Huang, et al. Apneumonia outbreak associated with a new coronavirus of probable bat origin.
Nature , 579(7798):270–273, 2020.[2] J. K. Y. Yap, M. Moriyama, and A. Iwasaki. Inflammasomes and pyroptosis as therapeutic targets forCOVID-19.
J. Immunol. , Jun 2020.[3] M. Z. Tay, C. M. Poh, L. R´enia, P. A. MacAry, and L. F. P. Ng. The trinity of COVID-19: Immunity,inflammation and intervention.
Nat. Rev. Immunol. , 20(6):363–374, 06 2020.[4] M. Soy, G. Keser, P. Atag¨und¨uz, F. Tabak, I. Atag¨und¨uz, and S. Kayhan. Cytokine storm in COVID-19:Pathogenesis and overview of anti-inflammatory agents used in treatment.
Clin. Rheumatol. , 39(7):2085–2094, Jul 2020.[5] D. Tang, P. Comish, and R. Kang. The hallmarks of COVID-19 disease.
PLoS Pathog. , 16(5):e1008536,05 2020.[6] A. Shah. Novel coronavirus-induced NLRP3 inflammasome activation: A potential drug target in thetreatment of COVID-19.
Front. Immunol. , 11, 2020.[7] M. Z. Ratajczak and M. Kucia. SARS-CoV-2 infection and overactivation of NLRP3 inflammasome asa trigger of cytokine “storm” and risk factor for damage of hematopoietic stem cells.
Leukemia , pages1–4, 2020.[8] Y. Fu, Y. Cheng, and Y. Wu. Understanding SARS-CoV-2-mediated inflammatory responses: Frommechanisms to potential therapeutic tools.
Virol. Sin. , pages 1–6, 2020.[9] S. Nagashima, M. C. Mendes, A. P. C. Martins, N. H. Borges, T. M. Godoy, A. F. M. Ribeiro, F. daSilva Deziderio, L. de Noronha, and C. Machado-Souza. The endothelial dysfunction and pyroptosisdriving the SARS-CoV-2 immune-thrombosis.
MedRxiv , 2020.[10] Y. Jamilloux, T. Henry, A. Belot, S. Viel, M. Fauter, T. El Jammal, T. Walzer, B. Fran¸cois, andP. S`eve. Should we stimulate or suppress immune responses in COVID-19? Cytokine and anti-cytokineinterventions.
Autoimmun. Rev. , 19(7):102567, Jul 2020.[11] C. Y. Taabazuing, M. C. Okondo, and D. A. Bachovchin. Pyroptosis and apoptosis pathways engage inbidirectional crosstalk in monocytes and macrophages.
Cell Chem. Biol. , 24(4):507–514, Apr 2017.[12] N. Kelley, D. Jeltema, Y. Duan, and Y. He. The NLRP3 inflammasome: an overview of mechanisms ofactivation and regulation.
Int. J. Mol. Sci. , 20(13):3328, 2019.[13] P. Song, W. Li, J. Xie, Y. Hou, and C. You. Cytokine storm induced by SARS-CoV-2.
Clin. Chim.Acta , 509:280–287, Jun 2020.[14] C. Huang, Y. Wang, X. Li, L. Ren, J. Zhao, Y. Hu, L. Zhang, G. Fan, J. Xu, X. Gu, et al. Clinicalfeatures of patients infected with 2019 novel coronavirus in Wuhan, China.
The Lancet , 395(10223):497–506, 2020.[15] M. Bertinaria, S. Gastaldi, E. Marini, and M. Giorgis. Development of covalent NLRP3 inflammasomeinhibitors: Chemistry and biological activity.
Arch. Biochem. Biophys. , 670:116–139, 07 2019.[16] K. Tsuchiya. Inflammasome-associated cell death: Pyroptosis, apoptosis, and physiological implications.
Microbiol. Immunol. , 64(4):252–269, Apr 2020.[17] S. L. Fink and B. T. Cookson. Apoptosis, pyroptosis, and necrosis: Mechanistic description of deadand dying eukaryotic cells.
Infect. Immun. , 73(4):1907–1916, Apr 2005.August 11, 2020 12/2218] M.P. Lythgoe and P. Middleton. Ongoing clinical trials for the management of the COVID-19 pandemic.
Trends Pharmacol Sci.
Mater Meth-ods , 10:2869, 2020.[20] A. Stutz, D. T. Golenbock, and E. Latz. Inflammasomes: Too big to miss.
J. Clin. Invest. , 119(12):3502–3511, Dec 2009.[21] S. Christgen, D. E. Place, and T. D. Kanneganti. Toward targeting inflammasomes: Insights into theirregulation and activation.
Cell Res. , 30(4):315–327, Apr 2020.[22] Z. Zheng and G. Li. Mechanisms and therapeutic regulation of pyroptosis in inflammatory diseases andcancer.
Int. J. Mol. Sci. , 21(4), Feb 2020.[23] M. G. Sanders, M. J. Parsons, A. G. Howard, J. Liu, S. R. Fassio, J. A. Martinez, and L. Bouchier-Hayes. Single-cell imaging of inflammatory caspase dimerization reveals differential recruitment toinflammasomes.
Cell Death Dis , 6:e1813, Jul 2015.[24] C. Semino, S. Carta, M. Gattorno, R. Sitia, and A. Rubartelli. Progressive waves of IL-1 β release byprimary human monocytes via sequential activation of vesicular and gasdermin D-mediated secretorypathways. Cell Death Dis. , 9(11):1–14, 2018.[25] G. Lopez-Castejon and D. Brough. Understanding the mechanism of IL-1 β secretion. Cytokine GrowthFactor Rev. , 22(4):189–195, 2011.[26] P. Broz, P. Pelegr´ın, and F. Shao. The gasdermins, a protein family executing cell death and inflam-mation.
Nat. Rev. Immunol. , pages 1–15, 2019.[27] L. DiPeso, D. X. Ji, R. E. Vance, and J. V. Price. Cell death and cell lysis are separable events duringpyroptosis.
Cell Death Disc. , 3(1):1–10, 2017.[28] D. Brough and N. J. Rothwell. Caspase-1-dependent processing of pro-interleukin-1 β is cytosolic andprecedes cell death. J. Cell Sci. , 120(5):772–781, 2007.[29] Y. Huang, H. Jiang, and et al. Chen, Y. Tranilast directly targets NLRP3 to treat inflammasome-drivendiseases.
EMBO Mol. Med. , 10(4):e8689, 2018.[30] K. Schleich and I. N. Lavrik. Mathematical modeling of apoptosis.
Cell Comm. Signal. , 11(1):1–7, 2013.[31] S. L. Spencer and P. K. Sorger. Measuring and modeling apoptosis in single cells.
Cell , 144(6):926–939,2011.[32] W. Wang and T. Zhang. Caspase-1-mediated pyroptosis of the predominance for driving CD4 T cellsdeath: A nonlocal spatial mathematical model.
Bull. Math. Biol. , 80(3):540–582, 2018.[33] D. Veltman, T. Laeremans, E. Passante, and H. J. Huber. Signal transduction analysis of the NLRP3-inflammasome pathway after cellular damage and its paracrine regulation.
J. Theor. Biol. , 415:125–136,2017.[34] Y. Bozkurt, A. Demir, B. Erman, and A. G¨ul. Unified modeling of familial mediterranean fever andcryopyrin associated periodic syndromes.
Comp. Math. Meth. Med. , 2015, 2015.[35] F. L´opez-Caamal and H. J. Huber. Stable IL-1 β -activation in an inflammasome signalling model dependson positive and negative feedbacks and tight regulation of protein production. IEEE ACM T. Comput.Bi. , 16(2):627–637, 2018.[36] A. V. Bagaev, A. Y. Garaeva, E. S. Lebedeva, A. V. Pichugin, R. I. Ataullakhanov, and F. I.Ataullakhanov. Elevated pre-activation basal level of nuclear NF- κ B in native macrophages acceler-ates LPS-induced translocation of cytosolic NF- κ B into the cell nucleus.
Sci. Reports , 9(1):1–16, 2019.August 11, 2020 13/2237] M. S. Salahudeen and P. S. Nishtala. An overview of pharmacodynamic modelling, ligand-bindingapproach and its application in clinical practice.
Saudi Pharm. J. , 25(2):165–175, Feb 2017.[38] N. M. de Vasconcelos, N. Van Opdenbosch, H. Van Gorp, E. Parthoens, and M. Lamkanfi. Single-cellanalysis of pyroptosis dynamics reveals conserved GSDMD-mediated subcellular events that precedeplasma membrane rupture.
Cell Death Diff. , 26(1):146–161, 01 2019.[39] MATLAB. version 1.8.0 202 (R2019n) . The MathWorks Inc., Natick, Massachusetts, 2019.[40] K. Tsuchiya. Inflammasome-associated cell death: Pyroptosis, apoptosis, and physiological implications.
Microb. Immunol. , 64(4):252–269, 2020.[41] . He, H. Wan, L. Hu, P. Chen, X. Wang, Z. Huang, Z. Yang, C. Zhong, and J. Han. Gasdermin D is anexecutor of pyroptosis and required for interleukin-1 β secretion. Cell Res. , 25(12):1285–1298, 2015.[42] J. Rex, A. Lutz, L. E. Faletti, U. Albrecht, M. Thomas, J. G. Bode, C. Borner, O. Sawodny, andI. Merfort. IL-1 and TNF differentially influence NF- κ B activity and FasL-induced apoptosis in primarymurine hepatocytes during LPS-induced inflammation.
Front. Physiol , 10:117, 2019.[43] Y. Wang, G. An, A. Becker, C. Cockrell, N. Collier, M. Craig, C. L. Davis, J. Faeder, A. N. F. Versypt,J. F. Gianlupi, et al. Rapid community-driven development of a SARS-CoV-2 tissue simulator.
BioRxiv ,2020.[44] A. Ghaffarizadeh, R. Heiland, S. H. Friedman, S. M. Mumenthaler, and P. Macklin. PhysiCell: An opensource physics-based cell simulator for 3-D multicellular systems.
PLoS Comput. Biol. , 14(2):e1005991,02 2018.[45] J. S. Lolkema and D. Slotboom. The Hill analysis and co-ion–driven transporter kinetics.
J. Gen.Physio. , 145(6):565–574, 2015.[46] C. Juliana, T. Fernandes-Alnemri, S. Kang, A. Farias, F. Qin, and E. S. Alnemri. Non-transcriptionalpriming and deubiquitination regulate NLRP3 inflammasome activation.
J. Biol. Chem. , 287(43):36617–36622, 2012.[47] S. Han, T. B. Lear, J. A. Jerome, S. Rajbhandari, C. A. Snavely, D. L. Gulick, K. F. Gibson, C. Zou,B. B. Chen, and R. K. Mallampalli. Lipopolysaccharide primes the NALP3 inflammasome by inhibitingits ubiquitination and degradation mediated by the SCFFBXL2 E3 ligase.
J. Bio. Chem. , 290(29):18124–18133, 2015.[48] M. A. Moors and S. B. Mizel. Proteasome-mediated regulation of interleukin-1 β turnover and exportin human monocytes. J. Leukocyte Biol. , 68(1):131–136, 2000.[49] F. Mart´ın-S´anchez, C. Diamond, M. Zeitler, A. I. Gomez, A. Baroja-Mazo, J. Bagnall, D. Spiller,M. White, M. J. D. Daniels, A. Mortellaro, et al. Inflammasome-dependent IL-1 β release depends uponmembrane permeabilisation. Cell Death Diff. , 23(7):1219–1231, 2016.August 11, 2020 14/22
Appendix
A.1 Complete system of ordinary differential equations
Here, we provide the full system of ODEs, as described in Sections 2.1 through to 2.8. Terms written inburgundy pertain to drug dynamics and can be excluded from the model if one is interested in simulatingthe pyroptosis pathway in the absence of drugs.[NF- κ B c ] dt = − S F k a [NF- κ B c ] + k b [NF- κ B n ] , [NF- κ B n ] dt = S F k a [NF- κ B c ] − k b [NF- κ B n ] ,d [NLRP3 i ] dt = α [NF- κ B n ] γ NF NF γ NF + [NF- κ B n ] γ NF − S k [NLRP3 i ] − δ [NLRP3 i ] ,d [NLRP3 a ] dt = S k [NLRP3 i ] − k F [NLRP3 a ] − δ [NLRP3 a ] − k +TR [TR] [NLRP3 a ] + k − TR [TR · NLRP3 a ] ,d [NLRP3 o ] dt = k F [NLRP3 a ] ,d [ASC f ] dt = − k (cid:0) − F (cid:1) [NLRP3 o ] [ASC f ] ,d [ASC b ] dt = k (cid:0) − F (cid:1) [NLRP3 o ] [ASC f ] ,d [pro-C1] dt = − k [ASC b ] [pro-C1] ,d [ C dt = k [ASC b ] [pro-C1] ,d [GSDMD] dt = − α [C1] γ C C1 γ C + [C1] γ C [GSDMD] ,d [GSDMD-N] dt = α [C1] γ C C1 γ C + [C1] γ C [GSDMD] ,d [pro-IL-1 β ] dt = α [NF- κ B n ] γ NF NF γ NF + [NF- κ B n ] γ NF − α [C1] γ C C1 γ C + [C1] γ C [pro-IL-1 β ] − δ [pro-IL-1 β ] ,d [IL-1 β c ] dt = α [C1] γ C C1 γ C + [C1] γ C [pro-IL-1 β ] − k G [IL-1 β c ] − δ [IL-1 β c ] ,d [IL-1 β e ] dt = k G [IL-1 β c ] ,d [pro-IL-18] dt = − α [C1] γ C C1 γ C + [C1] γ C [pro-IL-18] ,d [IL-18 c ] dt = α [C1] γ C C1 γ C + [C1] γ C [pro-IL-18] − k G [IL-18 c ] ,d [IL-18 e ] dt = k G [IL-18 c ] ,dVdt = k G V,d [TR] dt = − k +TR [TR] [NLRP3 a ] + k − TR [TR · NLRP3 a ] ,d [TR · NLRP3 a ] dt = k +TR [TR] [NLRP3 a ] − k − TR [TR · NLRP3 a ] . Where, S , F , S and G are respectively defined in Equations (4), (5), (9) and (22).August 11, 2020 15/22 .2 Abbreviations used in text • ARDS: acute respiratory distress syndrome • ASC: apoptosis-associated speck-like protein containing a CARD • CARD: caspase activation and recruitment domain • COVID-19: coronavirus disease 2019 • CD: Catalytic domain • C1-CD: the catalytic domain of caspase-1 • DAMP: damage associated molecular pattern • GSDMD: gasdermin D • GSDMD-N: N-terminal domain of GSDMD • IL-1 β : interleukin 1 β • IL-18: interleukin 18 • LLR: leucine-rich repeat • NBD: nucleotide-binding and oligomerisation domain • NF- κ B: nuclear factor kappa-light-chain-enhancer of activated B cells • NLRP3: NBD, LRR-containing receptors with an N-terminal PYD 3 • ODE: ordinary differential equation • PAMP: pathogen associated molecular pattern • pro-IL-1 β : pro-form of interleukin 1 β • pro-IL-18: pro-form of interleukin 18 • PYD: pyrin domain • SARS-CoV-2: severe acute respiratory syndrome coronavirus 2 • TLR: toll like receptor
A.3 Detailed model description and motivation
In this section we provide a detailed model description and justification. Model variables/components arelisted in Table 1 and model parameters are discussed in Appendix A.4.
Signal 1:
The initial signal, S , as defined in equation (4), is here considered to be either on , S = 1, or off , S = 0. This binary choice is motivated by the modelling assumption that there will be no inslammasome-inducing signalling until the DAMP/PAMP-sensing TLR activity is high enough. Therefore any signal levelbelow some threshold value, here e.g., S = 1, would not stimulate a downstream response. We here alsoassumed that once S is turned on , this signal is kept at a constant level until the cell dies. In future modeladaptations this binary and constant signalling can be modified to allow for continuous responses determinedby internal and external stimuli.August 11, 2020 16/22ymbol Component Description[ASC b ]( t ) concentration of bound ASC protein[ASC f ]( t ) concentration of free cytosolic ASC protein[C1]( t ) concentration of activated caspase-1[GSDMD]( t ) concentration of inactive GSDMD[GSDMD-N]( t ) concentration of activated GSDMD-N[IL-1 β c ]( t ) concentration of active IL-1 β in the cytoplasm[IL-1 β e ]( t ) concentration of active IL-1 β external[IL-18 c ]( t ) concentration of active IL-18 in the cytoplasm[ IL − e ]( t ) concentration of active IL-18 external[NF- κ B c ]( t ) concentration of cytosolic NF- κ B[NF- κ B n ]( t ) concentration of nuclear NF- κ B[NLRP3 a ]( t ) concentration of active NLRP3 protein[NLRP3 i ]( t ) concentration of inactive NLRP3 protein[NLRP3 o ]( t ) concentration of oligomerised NLRP3 protein[pro-C1]( t ) concentration of pro-caspase-1[pro-IL-1 β ]( t ) concentration of pro-IL-1 β [pro-IL-18]( t ) concentration of pro-IL-18 S first priming signal (e.g., virus) S second activation signal (e.g., ATP)[TR]( t ) concentration of free tranilast (drug)[TR · NLRP3 a ]( t ) concentration of tranilast - NLRPR a complex V ( t ) cell volume Table 1.
The symbols and descriptions of the components considered in the model, in alphabetical order.The bracket notation [ ] denotes compound concentration. Implicit time-dependence is denoted by (t).
NF- κ B dynamics:
The dynamics of the protein concentrations of cytoplasmic NF- κ B, [NF- κ B c ], andnuclear NF- κ B, [NF- κ B n ], are described in equations (1) and (2), respectively. When S = 1, we assume areaction of the form, NF- κ B c k a (cid:10) k b NF- κ B n . (32)Applying the law of mass action, we obtain the ODEs for both concentrations [NF- κ B c ] and [NF- κ B n ]. NLRP3 dynamics:
The dynamics of the concentrations of inactive NLRP3, [NLRP3 i ], active NLRP3,[NLRP3 a ] and oligomerised NLRP3, [NLRP3 o ] are described in Equations (6), (7) and (8), respectively. Todescribe the transcription of inactive NLRP3, which is induced by [NF- κ B n ] levels, we use a Hill function [37],[NF- κ B n ] γ NF N F γ NF + [NF- κ B n ] γ NF . (33)Here, NF represents the level of [NF- κ B n ] required for the Hill function to reach the value of a half, oftenreferred to the half-max value of [NF- κ B n ] [45]. To incorporate activation of NLRP3 in our model we considerthe reaction, NLRP3 i k −→ NLRP3 a . (34)Furthermore, for both inactive and active NLRP3 we include decay reactions of the form,NLRP3 i δ −→ ∅ , NLRP3 a δ −→ ∅ , (35)where ∅ denotes material outwith the system modelled. Finally, to include oligomerisation of NLRP3, if F = 1, we consider the reactionNLRP3 a + NLRP3 a k −→ NLRP3 o + NLRP3 o . (36)August 11, 2020 17/22or all above reactions, we apply the law of mass action to obtain the ODEs for [NLRP3 i ], [NLRP3 a ] and[NLRP3 o ]. The NLRP3 a -to-NLRP3 o oligomerisation step is switched off once [NLRP3 o ] reaches the level n ,where n corresponds to the concentration required to form an inflammasome base. We include this switchsince only one inflammasome is formed per cell during the pyroptosis process [20]. In the model, this switch iscontrolled by the step function F which, as described in Equation (5), is equal to one prior to inflammasomebase formation and zero afterwards. Signal 2:
In the model we assume that the second signal, S , is either on or off . For all shown simulationresults, we assume that S is on , so that S = 1, for all time points in the simulation. Through experiments,Juliana et al. [46] showed that increasing the time that a cell was exposed to ATP (yielding a second signal S to turn on ) did not alter the subsequent inflammasome dynamics or the expression of involved proteins.Thus, it could be argued that assuming that S is constant, i.e., on or off , throughout the simulation is areasonable model assumption. ASC dynamics:
In the model, we that assume free ASC, ASC f , binds to the inflammasome base, i.e., oligomerised [NLRP3 o ] when F = 0, through a binding reaction,ASC f + NLRP3 o k −→ ASC b + NLRP3 o . (37)Again, mass action kinetics allows us to formulate this interaction in terms of ODEs as provided in Equations(10) and (11). Caspase-1 dynamics:
We assume that pro-caspase-1, pro-C1, binds to inflammasome-bound ASC, ASC b .This binding is enabled in the model once [ASC b ] levels are larger than zero through a the reaction,ASC b + pro-C1 k −→ ASC b + C1 . (38)Using mass action kinetics, this reaction can be formulated in terms of the ODEs in Equations (13) and (14). Gasdermin dynamics:
The dynamics of GSDMD and GSDMD-N are described in Equations (16) and(17). We use a [C1]-dependent Hill function to capture the fact that the cleavage of GSDMD is mediated bycapsase-1. From here, the Hill function inclusive ODE formulation follows the methodology previously de-scribed for NF- κ B-dependent transcription of inactive NLRP3, or in the ODEs, [NF- κ B]-dependent elevationof [NLRP3 i ]. Cytokine dynamics (IL-1 β ): The dynamics of IL-1 β in the pro-form, pro-IL-1 β , the active form inthe cytoplasm, IL-1 β c , and the active form released into the external environment, IL-1 β e are described inEquations (19), (20) and (21). We choose a [NF- κ B]-dependent Hill function to describe the transcriptionof pro-IL-1 β , as is done for the transcription of NLRP3. Further, to incorporate the cleavage of pro-IL-1 β by C1, we use a [C1]-dependent Hill function. That is, we assume similar binding-and-cleavage kinetics forGSDMD and the cytokines. Furthermore, for both concentrations [pro-IL-1 β ] and [IL-1 β c ], we include decayreactions of the form, pro-IL-1 β δ −→ ∅ , IL-1 β c δ −→ ∅ , (39)where ∅ denotes material outwith the system modelled. Finally, we consider the secretion of cytoplasmicIL-1 β to the external environment to be dependent on the GSDMD-N derived membrane pores. That is, weassume a [GSDMD-N]-dependent reaction of the form,IL-1 β c k −→ IL-1 β e . (40)Again, mass action kinetics allows us to write interactions pertaining to IL-1 β dynamics in terms of ODEs,as provided in Equations (19), (20) and (21).August 11, 2020 18/22 ytokine dynamics (IL-18): The dynamics of IL-18 in the pro-form, pro-IL-18, active form in thecytoplasm, IL-18 c and the active form released to the external environment IL-18 c are described in Equations(23), (24) and (25). Caspase-1 induced cleavage of pro-IL-18 is modelled in the same way as the cleavageof pro-IL-1 β . Moreover, the secretion of cytoplasmic IL-18 to the external environment is modelled as forIL-1 β so that we assume a reaction of the form,IL-18 c k −→ IL-18 e . (41)Mass action kinetics allows us to formulate IL-18 related interaction as ODEs, as given in Equations (23),(24) and (25). Volume of the cell:
The dynamics of the volume of the cell, V , are given in Equation (27). It is assumedthat once pores on the cell membrane form, water can enter the cell and cause the volume of the cell toincrease at the rate k . Data from single-cell analysis by de Vasconcelos et al. [38] demonstrates pyroptosingcells increasing 50% in volume prior to the complete membrane rupture. This volume increase was reportedto be gradual (hence here approximated as linear) and started 10-15 minutes before membrane rupture. Modified drug targeting model:
The influence of the drug tranilast is included in a modified versionof the model which results in additional terms for the dynamics of NLRP3 a given in Equation (28). Alsoincluded in the drug model is the dynamics of the drug tranilast, TR, and the complex TR · NLRP3 a , as isin Equations (29) and (30). Here, we consider a binding reaction of the form,TR + NLRP3 a k +TR (cid:10) k -TR TR · NLRP3 a . (42)Mass action kinetics allows us to write this interaction in terms of ODEs. Tranislast is a covalent inhibitor[15], and thus we could an extra step in the above reaction. However, in lieu of relevant data, we arecurrently considering only one reaction step for TR binding to NLRP3 o in an effort to not include extramodel parameters. A.4 Model parameter values
The parameters used to perform numerical simulations of the model are provided in Table 2, unless statedotherwise. In the following paragraphs we described the reasoning behind the parameter values chosen.
NLRP3 dynamics:
In [47], it was found that NLRP3 had a half-life of approximately 6 hrs when stim-ulated with lipopolysaccharide (LPS), which translates to a decay rate of approximately 0.002 µ m min − .Therefore we choose, δ = 0 .
002 min − , as we assume measure concentrations as dimensionless, arbitrary units. Cytokine dynamics:
In [48], it was found that pro-IL-1 β had a half-life of approximately 2.5 hrs, whichtranslates to a decay rate of approximately 0.004 µ m min − . Therefore we choose, δ = 0 .
004 min − , as we assume measure concentrations as dimensionless, arbitrary units. The timeline of pyroptosis:
We use data from three experimental studies [36,38,49] to further estimatea timeline for the events regarded in our mathematical model.Bagaev et al. [36] studied NF- κ B kinetics in primary macrophages subjected to bacterial lipopolysaccharide(LPS). They reported that the nuclear NF- κ B concentration peaked at 10 minutes post LPS activation, afterAugust 11, 2020 19/22ymbol Description Value α transcription rate of NLRP3 0.025 (a.u.) min − α cleavage rate of GDSMD 0.08 min − α transcription rate of pro-IL-1 β − α cleavage rate of pro-IL-1 β − α cleavage rate of pro-IL-18 0.8 min − C1 half-max value of caspase-1 required for cleavage 0.3 a.u. δ decay rate of NLRP3 0.002 min − δ decay rate of pro-IL-1 β and IL-1 β − γ C1 Hill function coefficient of caspase-1 mediated cleavage 2 γ NF Hill function coefficient of NF- κ B mediated transcription 2 k a rate at which NF- κ B enters nucleus 0.3 min − k b rate at which NF- κ B leaves nucleus 0.03 min − k activation rate of NLRP3 0.07 min − k oligomerisation rate of NLRP3 0.07 min − k rate for NLRP3 - ASC interaction 0.02 (a.u.) − min − k rate for ASC - (pro-caspase-1) interactions 0.04 (a.u.) − min − k transport rate of IL-1 β out of cell 0.8 min − k transport rate of IL-18 out of cell 0.8 min − k rate at which cell volume increases 0.1 min − k +TR forward rate constant for the drug-NLRP3 a reaction 0.0035 (a.u.) − min − k − TR reverse rate constant for the drug-NLRP3 a reaction 0.000035 min − n relative level of NLRP3 o required for the inflammasome base 1 a.u.NF half-max value of NF- κ B required for transcription 0.3 a.u. V c critical volume at which cell ruptures 1.5 a.u. Table 2.
Parameters values used within numerical simulations of the model. The justification for thesevalues is given in Appendix A.4. Here ‘a.u.’ refers to arbitrary units of concentration or volume.August 11, 2020 20/22hich it decreased to a half-maximal level in a gradual manner over 100 minutes.In a single-cell analysis study, de Vasconcelos et al. [38] investigated subcellular key events that definepyroptosis. In their article, the authors provide a multitude of data pertaining to bone marrow-derivedmacrophages (BMDMs) that undergo caspase-1-dependent pyroptosis upon Bacillus anthracis lethal toxin(LeTx) stimulation. Some data and information from de Vasconcelos et al. [38] is listed below, • Cell volume was reported to increase gradually for approximately 13 minutes prior to membrane rup-ture. The cell volume increases by 50% prior to membrane rupture. (See Figure 6a,b in de Vasconcelos et al. [38]). • Supernatants from cell cultures were assayed for LDH activity. LDH is a cytosolic protein that isreleased by pyroptotic cells, and LDH activity is commonly used as a marker for plasma rupture inexperimental settings. Data showed that near-maximal supernatant LDH values were reached 120 min-utes post LeTx stimulation (see Figure 9c in de Vasconcelos et al. [38]). From this we approximate thetime between signal 1 turning on in the model, and complete membrane rupture, to be approximately2 hours. • The influx of Ca occurs before total membrane permeabilisation in pyroptosis (see Figure 7 in deVasconcelos et al. [38]). Further data suggests that changes in mitochondrial activity occur between(45+12=56) and (21+9=30) minutes prior to complete membrane rupture (see Figure 1 combined withFigure 3 and Supplementary Figure 3 in de Vasconcelos et al. [38]). If we assume that Ca influxand changes in mitochondrial activity commence between the priming and activation events, we canestimate that the NLRP3 inflammasome base is formed somewhere between 56 to 30 minutes priorto cell rupture. In the model we use the average value (43 minutes) as a time-point for when theinflammasome starts forming ( i.e. when the concentration of NLRP3 o reaches the threshold value n ).In [49], Mart´ın-S´anchez et al. measure the levels of pro-IL-1 β and external IL-1 β over time in mouse bonemarrow-derived macrophages (BMDMs) after inflammasome activation. Their results highlighted that therelease of IL-1 β coincided with membrane permeability (pores forming), and eventually all of the pro-IL-1 β present at the start of the experiments was cleaved and released from the cell, with approximately 90%released within 120 minutes.Using the above data from Bagaev et al. [36], de Vasconcelos et al. [38] and Mart´ın-S´anchez et al. [49], weconstruct the following timeline for the key events in our model: • S ) turns on . •
10 minutes (start time + 10 minutes). The nuclear NF- κ B concentration peaks. After this time point,the cytoplasmic-to-nuclear translocation of NF- κ B ceases. •
77 minutes (end time - 43 minutes). Enough NLRP3 has been transcribed, synthesised and activatedto form the inflammasome base. •
107 minutes (end time - 13 minutes). The GSDMD-N induced pores start forming, allowing the outfluxof IL-1 β and IL-18, as well as the influx of extracellular water causing the cell volume to start increasing. •
120 minutes (end time). The cell membrane completely ruptures. Around 90% of the pro-IL-1 β andpro-IL-18 will have be cleaved and subsequently released.We choose the values of α − α , C1 , γ C1 , γ NF , k a , k b , k − k , NF and V c to obtain dynamics whichfollow this timeline. The values are displayed in Table 2.August 11, 2020 21/22 rug (tranilast) dynamics: We are currently considering drug concentrations to be constant during thesimulated time span, effectively describing an in vitro scenario without drug wash out. This could be alteredto include drug dosages varying in time if appropriate. Immunoblot data presented by Huang et al. [29],pertaining to LPS-primed bone marrow derived macrophages (BMDMs) treated with TR, demonstrate dose-dependent (25-100 µ M) drug responses, where the maximal inhibition was reached at at 100 µ M of the drug.With this motivation, we investigate drug dosages in increments of 25 (here arbitrary concentration units)and adjust the k − T R /k + T R ratio to be “completely effective” at 100 a.u.. At this stage in the modeldevelopment, we interpret this “completely effectiveness” as the membrane rupture being delayed so that itdoes not occur within 300 simulated minutes. Accordingly, we set k − TR /k +TR = 0 . a.u. . A.5
MATLAB code and set up of numerical simulations
The mathematical model was implemented in
MATLAB using the ODE solver ode15s . We solve the systemof equations with the initial conditions (ICs),[NF- κ B c ](0) = 0 . , [NF- κ B n ](0) = 0 . , [ASC f ](0) = [pro-C1](0) = [GSDMD](0) = [pro-IL-18](0) = V (0) = 1 , while all other variables start out with a zero initial condition, to represent a cell in homeostasis. The ini-tial conditions pertaining to NF- κ B are motivated by previous experiments by Bagaev et al. [36], reportinginitial nuclear NF- κ B concentrations to be between 10% and 30% (of the total amount of cellular NF- κ B),depending on cell line. For other variables with a non-zero IC, we set the value 1 to be the initial concen-tration as we are currently working with arbitrary units of concentration/volume and considering relativecompound concentrations using the conservation laws provided throughout the manuscript. When includingdrug dynamics in the model, the drug IC is set to the desired drug dosage. The simulation results providedin Figures 3 and 4 are obtained by running the code with the ICs described above, along with the parametersgiven in Table 2 and the signals S and S turned on .The MATLAB code includes two events that allows for parameter updates in the simulation. The firstsuch event captures the formation of the inflammasome, so that when the concentration [NLRP3 o ] reaches n ,the required concentration for inflammasome formation, the variable F , as defined in Equation (5), switchesfrom 1 to 0. The second event is terminal and captures the complete membrane rupture of the cell, so thatwhen the cell volume V reaches a critical volume V c , the simulation stops.The MATLAB code, including user instructions, is available on our project GitHub page: https://github.com/Pyroptosishttps://github.com/Pyroptosis