key: cord-0155843-b7i8be5y authors: Hamis, Sara J; Macfarlane, Fiona R title: A single-cell mathematical model of SARS-CoV-2 induced pyroptosis and the anti-inflammatory response to the drug tranilast date: 2020-08-10 journal: nan DOI: nan sha: 8beaeca90fb0370d3698980abadc5665dc312aa5 doc_id: 155843 cord_uid: b7i8be5y Pyroptosis is an inflammatory mode of cell death that contributes to the cytokine storm associated with severe cases of coronavirus disease 2019 (COVID-19). Central to pyroptosis induced by severe acute 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-inflammatory drugs. One such drug, tranilast, is currently being evaluated as a COVID-19 treatment target in a clinical trial. In this study, we present a single-cell mathematical model that captures the formation of the NLRP3 inflammasome, pyroptotic cell death and drug-responses to tranilast. The model is formulated in terms of a system of ordinary differential equations (ODEs) that describe the dynamics of proteins involved in pyroptosis. 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). COVID-19 is a respiratory illness induced by the coronavirus strain SARS-CoV-2 [1] . Most people infected by the virus experience mild symptoms, however, in cases of severe disease, life-threatening symptoms can manifest [2] . These divergent disease trajectories have been, largely, attributed to differences in immune responses of infected hosts [3] . When cells register the presence of SARS-CoV-2 virions, a multitude of host-protective responses are triggered. Infected cells transmit signals that recruit immune cells, such as monocytes, macrophages and T-cells, to the site of infection. These immune cells act to eliminate the virus from the body, and thus they attack infected cells in which virions may replicate. Current research has shown that, upon active virion replication and release, SARS-CoV-2 can induce pyroptosis in both epithelial cells and immune cells [3] [4] [5] [6] [7] [8] [9] . Pyroptosis is an inflammatory and rapid mode of cell death that is characterised by the secretion of pro-inflammatory cytokines, cell swelling and, ultimately, membrane rupture resulting in the release of cytoplasmic contents to the extracellular environment [10, 11] . Furthermore, cytokines released by pyroptosing cells have been shown to induce pyroptosis in neighbouring bystander cells which can register the cytokines' damage associated molecular patterns (DAMPs) [12] . Fundamentally, both the recruitment of immune cells and pyroptotic cell death act to protect the host from the virus. Indeed, in hosts with healthy immune systems, virus-specific T-cells are recruited to the site of infection [3] , and cytokine levels are kept controlled by negative feedback regulations [13] . However, if the immune system contesting a viral infection is malfunctioning, a wide array of cytokines may be over-produced due to deregulation of the negative feedback that controls cytokine levels in healthy immune systems [13] . Such an over-compensating immune response may lead to an uncontrolled cytokine activity, commonly referred to as a cytokine storm [3] . Elevated levels of both pro-inflammatory and anti-inflammatory cytokines have been observed in severe cases of COVID-19 [2, 9, 13, 14] , commonly manifesting with severe symptoms such as pneumonia and acute respiratory distress syndrome (ARDS), which may lead to multiple organ failure [13] . Hence, cytokine storms are associated with poor clinical COVID-19 outcomes [2, 4, 6] . Therefore, suppressing the onset of cytokine storms and the subsequent inflammation, without completely cancelling-out host-protective effects of the immune 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 of the immune system untouched. In fact, there exists a number of anti-inflammatory drugs that can inhibit pyroptosis and, through the inhibition of pyroptosis, alter the cell death mode from inflammatory to noninflammatory [15] . Non-inflammatory modes of cell death includes apoptosis, which is characterised by cell shrinkage with cell membrane integrity maintained throughout cell death, whereas pyroptosis is characterised by 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 tissue damage that it brings. Secondly, it has been shown that tissue factors released upon pyroptosis may initiate blood coagulation cascades, therefore inhibiting pyroptosis may reduce the risk of blood clotting in COVID-19 cases [2] . One drug that inhibits pyroptosis is tranilast, an anti-inflammatory drug that is traditionally used for treating inflammatory diseases, such as asthma, in Japan and South Korea [15] . Tranilast is currently being re-purposed for treating COVID-19 related symptoms in a clinical trial 1 [2, 18] . In this subsection we describe key aspects of the intracellular pathway that leads to pyroptosis. For a more comprehensive description of the mechanisms driving pyroptosis, we refer the interested reader to detailed reviews [15, [19] [20] [21] [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 mathematical model. Black arrows represent reactions that involve mass transfer between depicted compounds. Blue arrows represent facilitation of compound formation or reactions. Red arrows represent signals that can be turned on or off in the model. Membrane pores (shown in red) induced by GSDMD-N allows for the outflux of inflammatory cytokines and the influx of extracellular water, causing the cell to swell and the cell membrane to eventually rupture. Abbreviations are listed in Appendix A.2. Central to pyroptosis is the formation of the inflammasome, a multi-protein complex consisting of three molecular units (a sensor, an adaptor and an executor) that together enable the release of cytokines and the ultimate membrane rupture that characterise pyroptosis. It has been shown that the NLRP3 inflammasome can 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 with an N-terminal pyrin domain (PYD) 3. Thus in this study we focus our attention on the inflammasome comprising 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 recruitment domain (CARD) and a PYD domain, while caspase-1 contains a CARD and a catalytic domain (CD). As is illustrated in Figure 1 , the inflammasome forms via homotypic interactions between NLRP3, ASC and pro-caspase-1. More specifically, NLRP3-ASC and ASC-caspase-1 bindings are facilitated by PYD-PYD and CARD-CARD interactions, respectively. Inflammasome formation can be triggered by the presence of microbial pathogens in a cell, or danger signals from other infected cells. Note that, upon formation, only one inflammasome is formed per cell [20] . At homeostasis, cells do not contain enough NLRP3 to produce an inflammasome base [21] . Instead, NLRP3 levels are elevated in cells when toll like receptors (TLRs) sense DAMPs or pathogen associated molecular patterns (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 activated B cells (NF-κB) is translocated to the nucleus, initiating the transcription and, by extension, the synthesis of (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 altered calcium signalling, potassium efflux and the generation of reactive oxygen species (ROS) [15] . In turn, active NLRP3 molecules can oligomerise and bind together to form a wheel shaped structure, constituting the inflammasome base. Through homotypic binding, ASC molecules can then bind to the NLRP3 inflammasome base. Thereafter, pro-caspase-1 can bind to ASC, enabling the proximity-induced dimerisation, and thereby the activation, of caspase-1 [23] . Active caspase-1 then acts to cleave the interleukins pro-IL-1β and pro-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 domain gasdermin D (GSDMD-N) which can form pores on the cell membrane. These pores enable the outflux of the inflamamtory cytokines IL-1β and IL-18 (in their mature forms) from the cytoplasm to the external environment [24, 25] . The released cytokines can recruit immune cells to the site of infection and initiate pyroptosis 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 membrane eventually 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 pyroptosis at homeostasis, and thus these proteins are not upregulated by the TLR-stimulated cytoplasm-to-nucleus translocation of NF-κB [19, 26] . Once the NLRP3 inflammasome is activated, cell lysis can be delayed but not prevented [27, 28] . However, if the formation of the NLRP3 inflammasome is inhibited, a series of intracellular events resulting in non-inflammatory cell death will instead be initiated. There exists multiple drugs that act to inhibit inflammasome formation in order to prevent pyroptosis, a rigorous list of covalent drugs that target the NLRP3 inflammasome can be found in a review by Bertinara et al. [15] . In this study, we include the pharmacodynamical effects of one of these drugs, namely, tranilast. Tranilast inhibits the formation of the NLRP3 inflammasome by covalently binding to the NBD domain of NLRP3 molecules, thus preventing NLRP3-NLRP3 interactions, and the subsequent formation of the NLRP3 inflammasome [29] . 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 specific aspects of inflammasome formation [33] [34] [35] . Within this work, we explicitly model the intracellular events that drive SARS-CoV-2 induced pyroptosis, from TLRs detecting DAMPs/PAMPs through to the ultimate membrane 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 in experimental studies. Thus, in this mathematical/computational work, we study the time evolution of nuclear NF-κB, the NLRP3 inflammasome and its components, the pore-forming protein GSDMD, the inflammatory cytokines IL-1β and IL-18, and cell swelling. We formulate a system of ODEs describing the dynamics of the key proteins involved in NLRP3 inflammasome 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 expand the model to include pharmacodynamic effects of the anti-inflammatory drug tranilast. Figure 1 provides a 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 modelling steps and parameters can be found in Appendices A.3 and A.4, respectively. Numerical simulation results are provided in Section 3. When TLRs register DAMPs or PAMPs, cytoplasmic NF-κB is translocated to the nucleus in order to initiate the transcription of inactive NLRP3 and pro-IL-1β. In the case of SARS-CoV-2, for example, these DAMPs and PAMPs can be induced by intracellular virion replication and release, or by cytokines from neighbouring, infected cells. In the model, we consider the cytoplasm-to-nucleus translocation process of NF-κB to be initiated by a signal S 1 that is turned on in response to TLR-sensed signals. The concentrations of cytoplasmic and nuclear NF-κB are here denoted by [NF-κB c ] and [NF-κB n ], respectively, and the rate constants are k 1a and k 1b . Once the inflammasome base has formed we assume that the translocation of NF-κB into the nucleus is stopped. Therefore, once F = 0 as defined in Equation (5), the translocation into 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 where we consider the total concentration of NF-κB to be constant over time [36] so that The DAMP/PAMP initiated TLR-signal S 1 is modelled as a binary on/off function such that Upon DAMP/PAMP induced TLR-signalling and nuclear NF-κB translocation, inactive NLRP3, NLRP3 i , is transcribed and subsequently synthesised. As the transcription and synthesis of inactive NLRP3 is promoted by nuclear NF-κB, we describe the production rate of NLRP3 using a Hill function. Specifically, a Hill function with a constant coefficient rate α 1 , Hill coefficient γ NF and half-max concentration-value NF 50 [37] . Inactive NLRP3 can self-activate in response to a secondary signal, S 2 , which can be turned on by multiple 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 inflammasome base, modelled by the concentration of oligomerised NLRP3, i.e., NLRP3 o . In the model, we thus consider NLRP3 protein concentrations in three different forms: Inactive and active NLRP3 decays in the model at a rate δ 1 . The forward reactions from the inactive-to-active, and activeto-oligomerised NLRP3 forms are here assumed to be irreversible, and are respectively described using the rates k 2 and k 3 . We additionally make the assumption that the oligomerisation of active NLRP3 will cease August 11, 2020 5/22 once an adequate amount, n, of NLRP3 has oligomerised and formed a base for the inflammasome, thus we introduce the step-function F , defined as We now incorporate all the above mechanisms to describe the rate of change of where the binary signal S 2 is set as Once the inflammasome base is formed, that is when F = 1 as defined in Equation (5) Note that, free ASC levels are adequate for pyroptosis at homeostasis [15, 19] , and thus the total amount of ASC is constant, Pro-caspase-1 can bind to inflammasome-bound ASC upon availability and subsequently dimerise into its activated, mature form, caspase-1. In the model, concentrations of pro-caspase-1 and caspase-1 are denoted by [pro-C1] and [C1], respectively. The activation occurs at a rate k 5 and is here assumed to be irreversible so that the rate of change of [pro-C1] and [C1] can be described by 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 GSDMD-N, which is formed as a result of caspase-1 cleavage of GSDMD, produces the pores on the cell membrane that are central to the pyroptosis process. In the model, we consider the concentrations of both GSDMD and GSDMD-N, denoted as [GSDMD] and [GSDMD-N], respectively. We describe the caspase-1 facilitated cleavage using a Hill function, where γ C1 is the Hill coefficient, and C1 50 is the half-max [C1] value. We additionally assume a specific rate for [GSDMD] cleavage, α 2 . Therefore, we describe the rate of change of [GSDMD] and [GSDMD-N] in the model, through the equations, Adequate levels of GSDMD required for pyroptosis are available in the cell at homeostasis [26] , therefore we include a conservation law, 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 active caspase-1 is available, the pro-forms of the interleukins can be cleaved into their activated forms. Subsequently, once membrane pores have been formed in response to GSDMD-N activity, cytoplasmic interleukins escape to the extracellular region. We here consider the concentrations of IL-1β/18 in pro-, cytoplasmic and 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 α 4 and α 5 , respectively. Once cleaved, cellular IL-1β c and IL-18 c escape through the GSDMD-N derived pores at the rates k 6 and k 7 , respectively. The transcription and synthesis of pro-IL-1β is promoted by [NF-κB n ], therefore we describe this process using a Hill function with constant coefficient α 3 , Hill coefficient γ NF and the half-max concentration of NF-κB, NF 50 . We additionally consider the decay of IL-1β at the rate δ 2 . We incorporate the above mechanisms to describe the rate of change of [pro-IL-1β], [IL-1β c ] and [IL-1β e ] over time as, Here, G is a function that allows for [GSDMD-N]-dependent cytokine outflux, and water influx, such that, The levels of pro-IL-18 in the cell are kept at a homeostatic level [15, 19] The GSDMD-N derived membrane pores that allow for the outflux of mature interleukins, also allow for the influx 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 increases gradually [38] . Thus, in the model, we consider the volume of the cell, V , to increase at a rate k 8 once pores are formed on the cell membrane. We describe the rate of change of the cell volume by the equation, 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 the inflammasome formation by preventing NLRP3-NLRP3 interactions [29] . Thus in the model, we consider TR binding to NLRP3 a to form the complex TR · NLRP3 a . The forward and reverse rate constants for this reaction 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, The rate of change of [TR] and [TR · NLRP3 a ] can now be included in the model as, where the following conservation law holds, The model is implemented in MATLAB [39] , where the system of ODEs (1)-(31), given in full in Appendix A.1, are solved numerically using the function ode15s. A more detailed description of the numerical set-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. We first consider the model without any drug, that is, we simulate the system given by Equations (1) We next consider numerical simulations of the model including the effects of the drug tranilast, that is the 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 initial concentrations (dosages) of the drug. We display the time evolution of each model component in a separate subplot, where the graph colour corresponds to drug dosage. For all components, increasing levels of the drug push pyroptotic events further in time. Note that the results for the maximal tested drug dosage results in many of the processes driving pyroptosis not occurring in the time frame examined. The drug tranilast specifically inhibits the NLRP3 a -to-NLRP3 o oligomerisation process. Therefore, when tranilast is added to the system, we observe that as the drug concentration [TR] increases, so does the amount of time it 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 ultimate membrane rupture, are delayed. Recall that if pyroptosis is delayed in a dying cell, the cell death mode can be altered from inflammatory to non-inflammatory. Pyroptosis has been identified as a key mechanism involved in the cytokine storm and the elevated inflammation associated with severe cases of COVID-19 [2] . Consequently, pyroptosis has been suggested as a COVID-19 treatment target. In order to investigate the pathway of pyroptosis in further detail, we have in this study formulated a single-cell mathematical model that captures the key proteins and processes involved in this pro-inflammatory mode of cell death. Specifically, we model the process from DAMP/PAMP-induced TLR signalling to cell death and membrane rupture. The model is formulated in terms of a system of ODEs describing the dynamics of protein concentrations and cell volume. We further expand the model to include pharmacodynamic effects of the anti-inflammatory drug tranilast, and our numerical simulations demonstrate that increasing drug levels result in delayed NLRP3 inflammasome formation, and by extension, delayed cell rupture. Recently it has been found that if GSDMD-N levels in the cell are low, inflammasome activation can lead to apoptosis instead of pyroptosis as GSDMD-N inhibits caspase-3 action [11, 19, 22, 40, 41] . Therefore, by delaying the cleavage of GSDMD, the cell death mechanism may be switched from inflammatory to non-inflammatory. As we have shown in our simulations inclusion of the drug tranilast can delay the formation 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-specified cell subjected to a non-specified pyroptosis-inducing stimuli. Our model parameter values were motivated by data from a collection of previous studies with different experimental settings. To quantitatively capture pyroptosis in a specific setting, we would require access to comprehensive data to update parameter choices accordingly. Furthermore, comprehensive data would also allow us to evaluate absolute concentrations for each of the proteins included within the model. With experimentally validated parameter values, a meaningful sensitivity analysis could be performed. However, as described in the Acknowledgements, this work is part of the Rapid Assistance in Modelling the Pandemic (RAMP) initiative and thus we choose to make our current work available now. In regards to further applications of this work, the single cell model developed in this study could be extended to investigate other anti-inflammatory drugs that target the NLRP3 inflammasome, a comprehensive review of such drugs is provided by Bertinara et al. [15] . This single cell model could also be incorporated into multi-cell models in order to study the effect of pyroptosis on a cell population or tissue scale. This would be especially interesting in a multi-scale model incorporating both epithelial cells and immune cells, as the release of cytoplasmic cytokines such as IL-1β and IL-18 into the external environment has been linked to contribute to the cytokine storm. It has been shown that IL-1β can induce pyroptosis in neighbouring cells [12] and that IL-18 can induce interferon (IFN) γ production in neighbouring cells, which plays a role in the recruitment of immune cells [42] . Therefore, it would be interesting to study these bystander effects in a multi-cell setting. Accordingly, we are collaborating with an interdisciplinary coalition of mathematicians and biologists to model SARS-CoV-2 induced within-host dynamics [43] . This multi-scale model is implemented in the PhysiCell [44] tissue simulator, allowing for in silico investigations of COVID-19 related dynamics of epithelial cells, immune cells, virions and cytokines. 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 in burgundy pertain to drug dynamics and can be excluded from the model if one is interested in simulating the pyroptosis pathway in the absence of drugs. Where, S 1 , F , S 2 and G are respectively defined in Equations (4), (5), (9) and (22) . • ARDS: acute respiratory distress syndrome • ASC: apoptosis-associated speck-like protein containing a CARD • 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 In this section we provide a detailed model description and justification. Model variables/components are listed in Table 1 and model parameters are discussed in Appendix A.4. Signal 1: The initial signal, S 1 , as defined in equation (4), is here considered to be either on, S 1 = 1, or off, S 1 = 0. This binary choice is motivated by the modelling assumption that there will be no inslammasomeinducing signalling until the DAMP/PAMP-sensing TLR activity is high enough. Therefore any signal level below some threshold value, here e.g., S 1 = 1, would not stimulate a downstream response. We here also assumed that once S 1 is turned on, this signal is kept at a constant level until the cell dies. In future model adaptations this binary and constant signalling can be modified to allow for continuous responses determined by internal and external stimuli. (6), (7) and (8), respectively. To describe the transcription of inactive NLRP3, which is induced by [NF-κB n ] levels, we use a Hill function [37] , Here, NF 50 represents the level of [NF-κB n ] required for the Hill function to reach the value of a half, often referred to the half-max value of [NF-κB n ] [45] . To incorporate activation of NLRP3 in our model we consider the reaction, Furthermore, for both inactive and active NLRP3 we include decay reactions of the form, where ∅ denotes material outwith the system modelled. Finally, to include oligomerisation of NLRP3, if F = 1, we consider the reaction where n corresponds to the concentration required to form an inflammasome base. We include this switch since only one inflammasome is formed per cell during the pyroptosis process [20] . In the model, this switch is controlled by the step function F which, as described in Equation (5), is equal to one prior to inflammasome base formation and zero afterwards. Signal 2: In the model we assume that the second signal, S 2 , is either on or off. For all shown simulation results, we assume that S 2 is on, so that S 2 = 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 2 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 2 is constant, i.e., on or off, throughout the simulation is a reasonable 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, 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, 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 by capsase-1. From here, the Hill function inclusive ODE formulation follows the methodology previously described for NF-κB-dependent transcription of inactive NLRP3, or in the ODEs, [NF-κB]-dependent elevation of [NLRP3 i ]. Cytokine dynamics (IL-1β): The dynamics of IL-1β in the pro-form, pro-IL-1β, the active form in the cytoplasm, IL-1β c , and the active form released into the external environment, IL-1β e are described in Equations (19) , (20) and (21) where ∅ denotes material outwith the system modelled. Finally, we consider the secretion of cytoplasmic IL-1β to the external environment to be dependent on the GSDMD-N derived membrane pores. That is, we assume a [GSDMD-N]-dependent reaction of the form, 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). The dynamics of IL-18 in the pro-form, pro-IL-18, active form in the cytoplasm, 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 cleavage of pro-IL-1β. Moreover, the secretion of cytoplasmic IL-18 to the external environment is modelled as for IL-1β so that we assume a reaction of the form, 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 assumed that once pores on the cell membrane form, water can enter the cell and cause the volume of the cell to increase at the rate k 8 . Data from single-cell analysis by de Vasconcelos et al. [38] demonstrates pyroptosing cells increasing 50% in volume prior to the complete membrane rupture. This volume increase was reported to 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 version of the model which results in additional terms for the dynamics of NLRP3 a given in Equation (28) . Also included in the drug model is the dynamics of the drug tranilast, TR, and the complex TR · NLRP3 a , as is in Equations (29) and (30) . Here, we consider a binding reaction of the form, 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 are currently considering only one reaction step for TR binding to NLRP3 o in an effort to not include extra model parameters. The parameters used to perform numerical simulations of the model are provided in Table 2 , unless stated otherwise. 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 stimulated with lipopolysaccharide (LPS), which translates to a decay rate of approximately 0.002 µm min −1 . Therefore we choose, δ 1 = 0.002 min −1 , 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, which translates to a decay rate of approximately 0.004 µm min −1 . Therefore we choose, 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 estimate a 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, after which 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 define pyroptosis. In their article, the authors provide a multitude of data pertaining to bone marrow-derived macrophages (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 rupture. 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 is released by pyroptotic cells, and LDH activity is commonly used as a marker for plasma rupture in experimental settings. Data showed that near-maximal supernatant LDH values were reached 120 minutes post LeTx stimulation (see Figure 9c in de Vasconcelos et al. [38] ). From this we approximate the time between signal 1 turning on in the model, and complete membrane rupture, to be approximately 2 hours. • The influx of Ca 2+ occurs before total membrane permeabilisation in pyroptosis (see Figure 7 in de Vasconcelos 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 with Figure 3 and Supplementary Figure 3 in de Vasconcelos et al. [38] ). If we assume that Ca 2+ influx and changes in mitochondrial activity commence between the priming and activation events, we can estimate that the NLRP3 inflammasome base is formed somewhere between 56 to 30 minutes prior to cell rupture. In the model we use the average value (43 minutes) as a time-point for when the inflammasome starts forming (i.e. when the concentration of NLRP3 o reaches the threshold value n). In [49] , Martín-Sánchez et al. measure the levels of pro-IL-1β and external IL-1β over time in mouse bone marrow-derived macrophages (BMDMs) after inflammasome activation. Their results highlighted that the release 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ánchez et al. [49] , we construct the following timeline for the key events in our model: • 0 minutes (start time). Signal 1 (S 1 ) 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 activated to form the inflammasome base. • 107 minutes (end time -13 minutes). The GSDMD-N induced pores start forming, allowing the outflux of 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β and pro-IL-18 will have be cleaved and subsequently released. We choose the values of α 1 − α 5 , C1 50 , γ C1 , γ NF , k 1a , k 1b , k 2 − k 8 , NF 50 and V c to obtain dynamics which follow this timeline. The values are displayed in Table 2 . A pneumonia outbreak associated with a new coronavirus of probable bat origin Inflammasomes and pyroptosis as therapeutic targets for COVID-19 The trinity of COVID-19: Immunity, inflammation and intervention Cytokine storm in COVID-19: Pathogenesis and overview of anti-inflammatory agents used in treatment The hallmarks of COVID-19 disease Novel coronavirus-induced NLRP3 inflammasome activation: A potential drug target in the treatment of COVID-19 SARS-CoV-2 infection and overactivation of NLRP3 inflammasome as a trigger of cytokine "storm" and risk factor for damage of hematopoietic stem cells Understanding SARS-CoV-2-mediated inflammatory responses: From mechanisms to potential therapeutic tools The endothelial dysfunction and pyroptosis driving the SARS-CoV-2 immune-thrombosis. MedRxiv Should we stimulate or suppress immune responses in COVID-19? Cytokine and anti-cytokine interventions Pyroptosis and apoptosis pathways engage in bidirectional crosstalk in monocytes and macrophages The NLRP3 inflammasome: an overview of mechanisms of activation and regulation Cytokine storm induced by SARS-CoV-2 Clinical features of patients infected with 2019 novel coronavirus in Wuhan Development of covalent NLRP3 inflammasome inhibitors: Chemistry and biological activity Inflammasome-associated cell death: Pyroptosis, apoptosis, and physiological implications Apoptosis, pyroptosis, and necrosis: Mechanistic description of dead and dying eukaryotic cells Ongoing clinical trials for the management of the COVID-19 pandemic Inflammasomes: Too big to miss Toward targeting inflammasomes: Insights into their regulation and activation Mechanisms and therapeutic regulation of pyroptosis in inflammatory diseases and cancer Single-cell imaging of inflammatory caspase dimerization reveals differential recruitment to inflammasomes Progressive waves of IL-1β release by primary human monocytes via sequential activation of vesicular and gasdermin D-mediated secretory pathways Understanding the mechanism of IL-1β secretion The gasdermins, a protein family executing cell death and inflammation Cell death and cell lysis are separable events during pyroptosis Caspase-1-dependent processing of pro-interleukin-1β is cytosolic and precedes cell death Tranilast directly targets NLRP3 to treat inflammasome-driven diseases Mathematical modeling of apoptosis Measuring and modeling apoptosis in single cells Caspase-1-mediated pyroptosis of the predominance for driving CD4 T cells death: A nonlocal spatial mathematical model Signal transduction analysis of the NLRP3-inflammasome pathway after cellular damage and its paracrine regulation Unified modeling of familial mediterranean fever and cryopyrin associated periodic syndromes Stable IL-1β-activation in an inflammasome signalling model depends on positive and negative feedbacks and tight regulation of protein production Elevated pre-activation basal level of nuclear NF-κB in native macrophages accelerates LPS-induced translocation of cytosolic NF-κB into the cell nucleus An overview of pharmacodynamic modelling, ligand-binding approach and its application in clinical practice Single-cell analysis of pyroptosis dynamics reveals conserved GSDMD-mediated subcellular events that precede plasma membrane rupture Inflammasome-associated cell death: Pyroptosis, apoptosis, and physiological implications Gasdermin D is an executor of pyroptosis and required for interleukin-1β secretion IL-1 and TNF differentially influence NF-κB activity and FasL-induced apoptosis in primary murine hepatocytes during LPS-induced inflammation Rapid community-driven development of a SARS-CoV-2 tissue simulator PhysiCell: An open source physics-based cell simulator for 3-D multicellular systems The Hill analysis and co-ion-driven transporter kinetics Non-transcriptional priming and deubiquitination regulate NLRP3 inflammasome activation Lipopolysaccharide primes the NALP3 inflammasome by inhibiting its ubiquitination and degradation mediated by the SCFFBXL2 E3 ligase Proteasome-mediated regulation of interleukin-1β turnover and export in human monocytes Inflammasome-dependent IL-1β release depends upon membrane permeabilisation The authors would like to thank the SARS-CoV-2 Tissue Simulation Coalition for ongoing collaboration and feedback on model development. Value α 1 transcription rate of NLRP3 0.025 (a.u.) min −1 α 2 cleavage rate of GDSMD 0.08 min −1 α 3 transcription rate of pro-IL-1β 0.007 (a.u.) min −1 α 4 cleavage rate of pro-IL-1β 0. 8 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 model development, we interpret this "completely effectiveness" as the membrane rupture being delayed so that it does not occur within 300 simulated minutes. Accordingly, we set k −TR /k +TR = 0.01a.u.. The mathematical model was implemented in MATLAB using the ODE solver ode15s. We while all other variables start out with a zero initial condition, to represent a cell in homeostasis. The initial conditions pertaining to NF-κB are motivated by previous experiments by Bagaev et al. [36] , reporting initial 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 concentration as we are currently working with arbitrary units of concentration/volume and considering relative compound concentrations using the conservation laws provided throughout the manuscript. When including drug dynamics in the model, the drug IC is set to the desired drug dosage. The simulation results provided in Figures 3 and 4 are obtained by running the code with the ICs described above, along with the parameters given in Table 2 and the signals S 1 and S 2 turned on.The MATLAB code includes two events that allows for parameter updates in the simulation. The first such 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), switches from 1 to 0. The second event is terminal and captures the complete membrane rupture of the cell, so that when 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/Pyroptosis.