A high order numerical method for solving Caputo nonlinear fractional ordinary differential equations

1.   Introduction

1.1.   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 infected by the virus experience mild symptoms, however, in severe cases of the disease, life-threatening symptoms can manifest ^[2]. These divergent disease trajectories have been, largely, attributed to differences in the immune response 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 into the extracellular environment ^[10,11]. Furthermore, interleukin-1$\beta$ (IL-1$\beta$), a cytokine released by pyroptosing cells, has been shown to induce pyroptosis in neighbouring bystander cells ^[12]. The cytokine interleukin-18 (IL-18), which is also released by pyroptosing cells, can act as a recruiter of immune cells ^{[13,14,15,16,17,18]}. Fundamentally, both the recruitment of immune cells and pyroptotic cell death act to protect the host from the virus. In hosts with healthy immune systems, virus-specific T-cells are recruited to the site of infection ^[3], and cytokine levels are kept under control by negative-feedback regulations ^[19]. 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 ^[19]. Such an over-compensating immune response may lead to 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,19,20], commonly manifesting with severe symptoms such as pneumonia and acute respiratory distress syndrome (ARDS), which may lead to multiple organ failure ^[19]. 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 being explored to combat symptoms of COVID-19 ^[19].

One approach to suppress cytokine storms is to inhibit pyroptosis, whilst leaving other functionalities of the host immune response untouched. In fact, there exists a number of anti-inflammatory drugs that may inhibit pyroptosis and, through this inhibition of pyroptosis, alter the mode of cell death from inflammatory to non-inflammatory ^[21]. Non-inflammatory modes of cell death include 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 ^[22,23]. Inhibiting pyroptosis may provide two key clinical advantages. Firstly, this inhibition could suppress the onset and effects of cytokine storms, and the resulting increased inflammation and tissue damage that they generate. 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].

1.2.   The pathway to pyroptosis

In this subsection we describe the key aspects of the intracellular pathway that leads to pyroptosis. Using biological details from this well-established pathway, we formulate a mathematical model in Section 2 that describes pyroptotic cell death. The model components are illustrated in Figure 1. For a more comprehensive description of the mechanisms driving pyroptosis, we refer the interested reader to detailed reviews ^{[17,21,24,25,26]}.

The formation of the NLRP3 inflammasome is a fundamental step in the pyroptosis pathway. This multi-protein complex consists of three molecular units (a sensor, an adaptor and an executor) and enables the cleavage the executor molecule caspase-1, which regulates pore formation and the release of cytokines in pyroptosis. Note that, only one inflammasome is formed per cell ^[17]. The inflammasome is named after its sensor molecule NLRP3, which is short for nucleotide-binding and oligomerisation domain (NBD) leucine-rich repeat (LRR)-containing receptors with an N-terminal pyrin domain (PYD) 3. In this study we focus our attention on the inflammasome comprising the sensor molecule NLRP3, the adaptor molecule apoptosis-associated speck-like protein (ASC), and caspase-1.

Homeostatic cells do not contain enough NLRP3 to produce an inflammasome base ^[25]. Instead, NLRP3 levels are elevated in cells when toll like receptors (TLRs) on the cell surface sense damage associated molecular patterns (DAMPs) or pathogen associated molecular patterns (PAMPs). As SARS-CoV-2 is a positive-sense RNA virus, it can be detected by TLRs and induce NLRP3 inflammasome formation ^[2,6,7]. Upon TLRs detecting DAMPs or PAMPs, the transcription factor NF-$\kappa$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$\beta$ (pro-IL-1$\beta$) are also regulated by NF-$\kappa$B. The increased transcription of (inactive) NLRP3 and pro-IL-1$\beta$ is often referred to as the priming step of pyroptosis, which is succeeded by the activation step ^[17].

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) ^[21]. In turn, active NLRP3 molecules can oligomerise, i.e., bind together, to form a wheel shaped structure, constituting the inflammasome base as depicted in Figure 2(a). Through homotypic binding, ASC molecules can then bind to the NLRP3 inflammasome base as illustrated in Figure 2(b). Thereafter, pro-caspase-1 can bind to ASC within the inflammasome, enabling the proximity-induced dimerisation, and thereby the activation, of caspase-1 ^[27]. Active caspase-1 then mediates the cleavage of the pro-interleukins, pro-IL-1$\beta$ and pro-IL-18, into their respective mature forms, IL-1$\beta$ and IL-18. Caspase-1 also cleaves the protein gasdermin D (GSDMD), releasing the active N-terminal domain of gasdermin D (GSDMD-N) which can form pores on the cell membrane. These pores enable the outflux of the inflammatory cytokines IL-1$\beta$ and IL-18 (in their mature forms) from the cytoplasm to the external environment ^[28,29]. The released cytokines can recruit immune cells to the site of infection and initiate pyroptosis in neighbouring cells. Furthermore, the GDSMD-N derived membrane pores also allow for the influx of extracellular material, e.g., water, into the cell. This influx causes cells to swell until their plasma membrane eventually ruptures, releasing cellular material into the extracellular region.

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-$\kappa$B ^[24,30]. Furthermore, once the NLRP3 inflammasome is activated, cell lysis can be delayed but not prevented ^[31,32]. However, if the formation of the NLRP3 inflammasome is delayed, a series of intracellular events resulting in non-inflammatory cell death may 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. ^[21]. In this study, we include the pharmacodynamical effects of a generic anti-inflammatory drug which inhibits the formation of the NLRP3 inflammasome by covalently binding to NLRP3 molecules to prevent NLRP3-NLRP3 interactions, as shown in Figure 2(a).

1.3.   Mathematical modelling of pyroptosis

The use of mathematical models to describe apoptotic cell death has been well established and reviewed ^[33,34]. However, there are significantly fewer mathematical models that describe pyroptotic cell death. Previous modelling works include more implicit descriptions of pyroptosis ^[35] and descriptions of specific aspects of inflammasome formation ^[36,37,38]. 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 a generic anti-inflammatory drug. The goal of our model is to capture the key aspects of the pyroptosis pathway that are commonly analysed in experimental studies. Thus, in this mathematical/computational work, we study the temporal evolution of NF-$\kappa$B, the NLRP3 inflammasome and its components, the pore-forming protein GSDMD, the pro-inflammatory cytokines IL-1$\beta$ and IL-18, and the cell volume.

1.4.   Mathematical modelling of within host-dynamics of COVID-19

Since COVID-19 was announced as a global pandemic ^[39], many researchers have developed and used mathematical tools and models to study the SARS-CoV-2 virus and COVID-19. These models can be split into two main categories, which are, external dynamics models (e.g., models that consider the spread of SARS-CoV-2 from person to person) and within-host models (e.g., models that capture the dynamics of SARS-CoV-2 and SARS-CoV-2 induced responses within the body). Agent-based models (ABMs) have been widely utilised in the past to model within-host viral dynamics ^{[40,41,42,43]}, with focuses on influenza ^{[44,45,46,47,48]}, oncolytic virotherapy ^[49,50] and inflammation ^[51,52]. ABMs can be used to describe how viruses spread within the body and circumvent the human immune system. Accordingly, several research groups are currently using multiscale ABMs to describe the within-host dynamics in response to SARS-CoV-2 ^[53,54,55]. As we describe in Section 3, our model of single-cell pyroptotic cell death is incorporated within a multiscale SARS-CoV-2 tissue simulator, that has been developed in response to the COVID-19 pandemic by an international and interdisciplinary research coalition. For the full details of this tissue simulator, we refer the reader to the project preprint ^[55].

2.   Mathematical model of pyroptosis in a single-cell

We formulate a system of ordinary differential equations (ODEs) describing the dynamics of the key components of NLRP3 inflammasome formation and pyroptosis. Specifically, the model includes the dynamics of NF-$\kappa$B (Section 2.1), NLRP3 (Section 2.2), ASC (Section 2.3), caspase-1 (Section 2.4), GSDMD (Section 2.5), IL-1$\beta$ and IL-18 (Section 2.6) and cell volume. Figure 1 provides a pictorial overview of the modelled pathway, and variable names are listed in Table B.1. Throughout this work, the bracket notation [ ] denotes molecule concentrations. Further details concerning modelling choices, reaction terms and parameterisation can be found in Appendices B and C. In Section 2.8, we expand the model to include the pharmacodynamic effects of an anti-inflammatory drug targeting NLRP3. We include the full system of equations and describe the set-up of the numerical simulations in Section 2.9. The results of these numerical simulations are described in Section 2.10, and we further investigate the model using sensitivity analysis in Section 2.11.

2.1.   NF-$\kappa$B dynamics

When TLRs register DAMPs or PAMPs, cytoplasmic NF-$\kappa$B is translocated to the nucleus in order to initiate the transcription of inactive NLRP3 and pro-IL-1$\beta$. In the case of SARS-CoV-2, these DAMPs and PAMPs can be induced by intracellular virion replication and release, or by cytokines from neighbouring, infected cells. Dynamics of NF-$\kappa$B translocation between the cytoplasm and the nucleus vary between situations and strongly depend on cell type and external stimuli. In recent work, Bagaev et al. ^[56] studied TLR-4 derived activation of NF-$\kappa$B in bone marrow derived macrophages (BMDMs) subjected to bacterial lipopolysaccharide (LPS). The authors found a sharp peak in NF-$\kappa$B nuclear translocation kinetics 10 minutes post LPS activation, after which the nuclear NF-$\kappa$B signalling gradually decreased until plateau. BMDMs stimulated with LPS have been shown to undergo pyroptosis ^[57,58], and we therefore aim to formulate a mathematical model that achieves nuclear NF-$\kappa$B levels that, over time, are $(i)$ initially rapidly increasing, $(ii)$ mono-peaking and $(iii)$ right-skewed. In order to achieve $(i-iii)$, we here choose to phenomenologically model nuclear NF-$\kappa$B levels using the function,

where [NF-$\kappa$B_n](0) is the amount of the cell's NF-$\kappa$B that is located in the nucleus prior to TLR-signalling. The DAMP/PAMP initiated TLR-signal, $S_1$, has been explicitly included as a binary on/off function such that,

In Eq (2.1), $h$ is a constant that denotes the peak elevation of the [NF-$\kappa$B_n] value, and $\tau$ denotes the time point (in units of minutes) at which the peak occurs. The variable $s$ controls the skewness of the [NK-$\kappa$B_n] peak over time. Accordingly, the parameters $h$, $\tau$ and $s$ can be varied to adjust the shape of the [NF-$\kappa$B_n] peak, as is shown in the Supplementary Material (Supplementary Material, S1).

2.2.   NLRP3 dynamics

Following nuclear NF-$\kappa$B translocation, inactive NLRP3, NLRP3_i, is transcribed and subsequently synthesised. As the transcription and synthesis of inactive NLRP3 is promoted by nuclear NF-$\kappa$B, we describe the production rate of NLRP3 using a Hill function ^[59]. Specifically, a Hill function with a constant coefficient rate $\alpha_{1}$, Hill coefficient $\gamma_{\rm{NF}}$ and half-max concentration-value NF₅₀. Inactive NLRP3 can become activated 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: [NLRP3_i], [NLRP3_a] and [NLRP3_o]. Inactive and active NLRP3 decays in the model at a rate $\delta_{1}$, however we do not include the decay of NLRP3_o in the model. The reason behind this choice of omitting NLRP3_o decay is that, in experimental works, the functions of the NLRP3 inflammasome are evident until cell death in the context of pyroptosis ^[60,61,62]. Furthermore, it is unclear whether disassembly of the inflammasome is possible in any context ^[63,64]. 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 the rates $k_{1}$ and $k_{2}$.

We now incorporate all the above mechanisms to describe the rate of change of [NLRP3_i], [NLRP3_a] and [NLRP3_o] over time as,

where [$\overline{\rm{NF- \mathit{ κ} B_n}}$] denotes the deviation from baseline in [NF-$\kappa$B_n] so that,

The binary signal $S_2$ is set as,

2.3.   ASC dynamics

Once the inflammasome base is formed, free ASC in the cell can bind to the NLRP3 inflammasome base. Thus in the model, we consider the concentrations of ASC in free and bound form, denoted [ASC_f] and [ASC_b], respectively. We consider this binding process to be irreversible and to occur at a rate $k_{3}$, hence the rate of change of [ASC_f] and [ASC_b] are here described by,

where the function $F([\rm{NLRP3_o}])$ is included in order to model the absence or presence of the NLRP3 inflammasome base. In order to ensure that binding of the ASC occurs only once the inflammasome base has formed in our continuous model, we approximate this two-state system as a continuous sigmoid function ^[65] such that,

Here, values $a = 1$ (a.u.) and $b = 2$ (a.u.) can be used to obtain $F([\rm{NLRP3_o}]) \approx 0$ for $[\rm{NLRP3_o}] < 1$ and $F(\rm[\rm{NLRP3_o}]) \approx 1$ otherwise. The value of $c$ is set to be large (e.g., 1000) in order to create a steep curve in the sigmoid. These values result in the binding of ASC once $[\rm{NLRP3_o}]$ reaches a threshold value $n\approx 1$ a.u., that is when oligomerised NLRP3 reaches this level the inflammasome is formed. {Note that the total amount of ASC is conserved in our model, }

Free ASC levels are adequate for pyroptosis at homeostasis ^[21,24]. Therefore, in an effort to minimise model complexity, we do not include the synthesis or decay of ASC in the model.

2.4.   Caspase-1 dynamics

Pro-caspase-1 can bind to inflammasome-bound ASC upon availability, and subsequently dimerise into its activated, mature form, caspase-1. In our model, concentrations of pro-caspase-1 and caspase-1 are denoted by [pro-C1] and [C1], respectively. Caspase-1 activation occurs at a rate $k_{4}$ and is here assumed to be irreversible so that the rate of change of [pro-C1] and [C1] can be described by,

so that,

Pro-caspase-1 is present in adequate levels in the cell prior to cellular responses to DAMPs or PAMPs ^[21,24]. Therefore, to minimise model complexity, neither the synthesis nor the decay of caspase-1 is included in the model.

2.5.   Gasdermin dynamics

GSDMD-N, which is formed as a result of the caspase-1 cleavage of GSDMD, produces 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 $\gamma_{\rm{C1}}$ is the Hill coefficient, and C1₅₀ is the half-max [C1] value. We additionally assume a specific rate for [GSDMD] cleavage, $\alpha_{2}$. Therefore, we describe the rate of change of [GSDMD] and [GSDMD-N] in the model, through the equations,

so that the following conservation law holds,

Adequate levels of GSDMD required for pyroptosis are available in the cell at homeostasis ^[30], therefore GSDMD synthesis and decay are omitted in the model.

2.6.   Cytokine dynamics

Translocation of NF-$\kappa$B, from the cytoplasm to the nucleus, induces the transcription and synthesis of pro-IL-1$\beta$, but does not up-regulate pro-IL-18 activity from homeostatic levels. 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 can be secreted into the extracellular region. We here consider the concentrations of IL-1$\beta$/18 in pro-, cytoplasmic and extracellular form, respectively denoted [pro-IL-1$\beta$/18], [IL-1$\beta$_c/18_c] and [IL-1$\beta$_e/18_e]. In the model, pro-IL-1$\beta$ and pro-IL-18 are cleaved by caspase-1 at rates $\alpha_{4}$ and $\alpha_{5}$, respectively. Once cleaved, cellular IL-1$\beta$_c and IL-18_c are secreted through the GSDMD-N derived pores at the rates $k_{5}$ and $k_{6}$, respectively.

IL-1$\beta$: The transcription and synthesis of pro-IL-1$\beta$ is promoted by [NF-$\kappa$B_n], in a similar way to NLRP3, therefore we describe this process using a Hill function with constant coefficient $\alpha_{3}$, Hill coefficient $\gamma_{\rm{NF}}$ and the half-max concentration of NF-$\kappa$B, NF₅₀. We additionally consider the decay of IL-1$\beta$ at the rate $\delta_{2}$. We incorporate the above mechanisms to describe the rate of change of [pro-IL-1$\beta$], [IL-1$\beta$_c] and [IL-1$\beta$_e] over time as,

where $[\overline{\rm{NF- \mathit{ κ} B_n}}]$ is defined in Eq (2.6) and $G$ is a non-dimensionalised value that allows for [GSDMD-N]-dependent cytokine outflux, and water influx, such that,

IL-18: We describe the rate of change of [pro-IL-18], [IL-18_c] and [IL-18_e] as,

so that,

The levels of pro-IL-18 in the cell are kept at a homeostatic level ^[21,24] and thus the synthesis and decay of IL-18 are not include in the model.

2.7.   Cell volume

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 its membrane eventually ruptures. Single-cell analysis has revealed that before the ultimate membrane rupture occurs, the cell volume increases gradually ^[60]. Thus, in the model, we consider the volume of the cell, $V$, to increase at a rate $k_{7}G$ 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}$, the cell ruptures and all cell processes cease.

2.8.   Modified drug targeting model

We consider a generic anti-inflammatory drug (Drug) which binds to active NLRP3 and thus inhibits the inflammasome formation by preventing NLRP3_a oligomerisation. To incorporate this in our model, we consider the drug binding to NLRP3_a to form the complex Drug $\cdot$ NLRP3_a. The forward and reverse rate constants for this reaction are denoted by $k_{+\rm{D}}$ and $k_{-\rm{D}}$, respectively. When including these drug mechanisms in our model, Equation (2.4) describing the dynamics of NLRP3_a must be modified to include two additional terms,

The rate of change of [Drug] and [Drug $\cdot$ NLRP3_a] can now be included in the model as,

where the following conservation law holds,

2.9.   Implementation of the pathway model

In order to reduce the pyroptosis model, as described in Sections 2.1-2.8, we impose the conservation laws (2.11), (2.14), (2.17), (2.25) and (2.30). We further set all total concentration constants to be 1 arbitrary unit of concentration (1 a.u.). Below, we provide the reduced system of equations,

with the following definitions,

The reduced model is implemented in MATLAB ^[66]}, where the system of ODEs, (2.31) with (2.32), is solved numerically using the built-in ODE-solver ode15s. Model outputs are measured in terms of arbitrary units of concentration or volume over time. The model parameters used in the simulations, and the justifications for these values, are provided in Table C.1 and Appendix C, respectively.

The initial conditions pertaining to NF-$\kappa$B are motivated by the works of Bagaev et al. ^[56], reporting initial nuclear NF-$\kappa$B concentrations to be between 10% and 30% (of the total amount of cellular NF-$\kappa$B), depending on cell line. Therefore we choose,

where the total amount of cellular NF-$\kappa$B is set to be 1 a.u.. Similarly, we set,

while all other initial concentrations are set to be 0 a.u.. For the results shown in Figure 3, we consider the model without drug intervention, whereas, in the results shown in Figure 4, we consider the initial concentration of the drug to be a value between 0 and 1 a.u.. Note, as we consider an in vitro situation, we consider the drug to be present from the start of the simulations and not decay within the simulation time.

A more detailed description of the numerical set-up can be found in Appendix D, which also includes instructions on how to access and run the code. Simulation results are provided and discussed in the next section.

2.10.   Results and discussion of the pathway model

We first consider the model without any drug, that is, we simulate the system of ODEs (2.31) with functions (2.32) and [Drug](0) = 0 a.u.. The simulation results are provided in Figure 3, in which the molecule concentrations (in arbitrary units) and cell volume are plotted over time. The results show that NF-$\kappa$B is rapidly translocated to the nucleus upon DAMP/PAMP-induced TLR signalling. After [NF-$\kappa$B_n] levels have peaked, the transcription, and by extension the synthesis, of NLRP3_i and pro-IL-1$\beta$ are reduced as NF-$\kappa$B leaves the nucleus. Synthesised NLRP3_i activates to become NLRP3_a, which in turn oligomerises and binds together to form the inflammasome base, here expressed in terms of the concentration of NLRP3_o. These dynamics are captured in the subplot displaying the time evolution of [NLRP3_i], [NLRP3_a] and [NLRP3_o]. Note that once the inflammasome base is formed, the oligomerised NLRP3 concentration plateaus and the binding of ASC to the inflammasome increases. This can be observed in the results, as [ASC_b] levels increase when [NLRP3_o] reaches the value $n$, where here $n = 1$ a.u.. Furthermore, as [ASC_b] levels increase, the dimerisation and activation of caspase-1 is facilitated so that [C1] levels increase. In turn, when [C1] levels increase, the cleavage of GSDMD, pro-IL-1$\beta$ and pro-IL-18 occurs. This is apparent in the subplots, as [GSDMD-N], [IL-1$\beta$_c] and [IL-18_c] levels start increasing once [C1] levels surpass zero a.u.. The outflux of mature interleukins and the influx of extracellular water to the cell are mediated by GSDMD-N derived pores. Thus, once GSDMD-N has been cleaved, so that [GSDMD-N] > 0 a.u., the cytoplasmic interleukin concentrations [IL-1$\beta$_c] and [IL-18_c] respectively decrease in favour of the extracellular concentrations [IL-1$\beta$_e] and [IL-18_e]. Increasing [GSDMD-N] levels also enable the cellular volume to increase, which eventually results in membrane rupture once the volume reaches the maximal capacity $V_c$, where here $V_c = 1.5$ a.u.. Note that when membrane rupture occurs, all intracellular mechanisms are ceased and thus the time progression in the subplots stops.

We next consider numerical simulations of the model including the effects of an anti-inflammatory drug, that is system of ODEs (2.31) with functions (2.32) and [Drug](0)> 0 a.u.. 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. These results show that increasing levels of the drug push pyroptotic events further forward 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 specifically inhibits the NLRP3 oligomerisation process. Therefore, when it is added to the system, we observe that as the drug concentration [Drug] 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.

2.11.   Sensitivity analyses of the pathway model

The mathematical model of the pyroptosis pathway includes the subset of model parameters $\Theta = \{ \theta_1 , \theta_2 , ..., \theta_p \}$, which are listed in Table C.1 in Appendix C. As is described in Appendix C, some model parameter values are directly obtained from published data ^[61,67], some are estimated from immunoblot data ^[58], whilst other values are collectively approximated through data fitting to time-course data that are available in the literature ^[56,60,68]. We note that such a data fitting approach may lead to large parameter uncertainties, and thus we conduct sensitivity analyses in order to investigate how sensitive the model outputs are to changes in the values of the model parameters $\Theta$, hereby referred to as the model inputs. Such analyses can help us understand how changes in the model input can affect the model output ^[69,70]. Therefore, this can allow us to identify which input parameters $\theta_j$ are the most influential, providing information on the importance of experimental retrieval of accurate parameters and possible model reductions. In this study, we investigate the system of ODEs (2.31) describing the pyroptosis pathway, using three different sensitivity analysis techniques, which are robustness analysis (Section 2.11.1), Latin hypercube analysis (Section 2.11.2) and a derivative-based method (Section 2.11.3). The molecule concentrations over time are the model outputs $y_i(t)$, where $i = 1, ..., q$, as labelled in Table B.1.

2.11.1.   Robustness analyses

Robustness analysis is a local sensitivity analysis technique in which one input parameter $\theta_j$ is varied, whilst the other $p-1$ parameters are kept fixed at their estimated values ^[69,70]. We compare the output molecule concentrations $y_i(t)$ over time for 7 different perturbed values of each parameter $\theta_j$, where the perturbed values range between $\pm20\%$ of the estimated value. In Figure 5, we include robustness analysis results for one input parameter, specifically $\alpha_1$, which denotes the rate coefficient for the transcription of inactive NLRP3. This figure shows that increasing the parameter value of $\alpha_1$ speeds up the inflammasome base formation and, thereby, the pyroptosis process. For investigated values of $\alpha_1$ below the threshold value of $\alpha_1 = 0.07$ (a.u.) min$^{-1}$, the pyroptosis process is delayed to an extent that NLRP3 formation does not occur within the simulated 500 minutes. For the estimated values of the parameter set $\Theta$, we consider such perturbations of $\alpha_1$ to be unfeasible, since we are here aiming to simulate a scenario in which the NLRP3 inflammasome base is formed at around 77 minutes, and the ultimate membrane rupture occurs at 120 minutes, as is described in Appendix C. Robustness analysis results for all model inputs $\theta_j, j = 1, ..., 21$ and model outputs $y_i(t)$, $i = 1, ..., 15$ are available in the Supplementary Material (Supplementary Material, S2). These analyses provides a visual overview of how increasing or decreasing each model parameter affects the output.

2.11.2.   Global sensitivity analysis

We perform a global sensitivity analysis, in which all parameters $\theta_j$ downstream of NF-$\kappa$B dynamics are simultaneously perturbed from their estimated values, using Latin hypercube sampling and analysis ^[69,70]. Results from the Latin hypercube analysis are provided in the form of scatter plots, where the model outputs are plotted over the investigated parameter ranges ($\pm 10\%$ of the estimated parameter values) for each input parameter. Here, the regarded output is the NLRP3_o concentration at 77 minutes (the time point where the inflammasome base is formed i.e., when [NLRP3_o] reaches the value $n = 1$ a.u. when using the estimated parameter values). The scatter plots in Figure 6 provide a visual means to study the relationship between three model inputs ($\alpha_1, k_1, \delta_1$) and the regarded output, for small perturbations of the input parameters. Scatter plots for the other input parameters are available in the Supplementary Material (Supplementary Material, S3). In order to quantify linear input-output associations, Figure 6 also includes Pearson correlation coefficients (R) for each input-output pair, evaluated over the regarded input parameter ranges ^[70]. A strong positive linear input-output correlation (R = 0.992) is found for input parameter $\alpha_1$, implying that, given the estimated parameter values of $\Theta$, our model results are especially sensitive to small perturbations of $\alpha_1$ which influences the transcription rate of inactive NLRP3. The linear input-output correlation is weakly positive for input parameter $k_1$ (R = 0.157), and negligible for $\delta_1$ (R = 0.011), within the investigated parameter ranges.

2.11.3.   Derivative-based sensitivity analysis

We perform a derivative-based sensitivity analysis, in which the equations describing the time evolution of the concentrations $y_i(t)$, where $i = 1, ..., q$, are differentiated with respect to each input parameter $\theta_j$ where $j = 1, ..., p$ ^[69,71]. We thus compute $q \times p$ partial derivatives $z^i_j(t) = \partial y_i(t) / \partial \theta_j$, at each time-point $t$ in which the system of ODEs (2.31) is numerically solved. These partial derivatives $z^i_j(t)$ are also referred to as the model sensitivities, and are here numerically computed using the Direct Method, as described by Dickinson and Gelinas ^[71]. The values $z^i_j(t)$ provide time-varying measurements of the sensitivity of the outputs $y_i(t)$ with respect to changes in values of the input parameters $\theta_j$. Results from the derivative-based sensitivity analysis are presented in Figure 7, in which output sensitivities with respect to $\alpha_1$, i.e, $z^i_{\alpha_1}$, are plotted over time. Results for sensitivities with respect to other model input parameters are available in the Supplementary Material (Supplementary Material, S4). Figure 7 shows that by increasing $\alpha_1$, the NLRP3_a levels increase at the start of the simulation (for times under $\sim$ 100 minutes). NLRP3_o levels also increase with $\alpha_1$, where the magnitude of this increase plateaus when the inflammasome base is formed. Thus the outputs [NLRP3_a] and [NLRP3_o] are both sensitive to variations of $\alpha_1$. As a downstream effect of this, the levels of C1 and GSDMD-N also increase with $\alpha_1$ which, in turn, alter the composition of pro-, cytoplasmic, and external IL-1$\beta$ over time, as well as alter the cell volume $V$.

3.   Integrating the pyroptosis pathway model with the PhysiCell framework

3.1.   Overview

Our single-cell pyroptosis model, described in Section 2, has been integrated within a multiscale tissue simulator that has been developed to investigate within-host dynamics in response to SARS-CoV-2 ^[55]. This work is part of a multidisciplinary, international coalition led by Paul Macklin (Indiana University, USA). The coalition is made up of several subetaoups who are developing individual submodels which are integrated to form a virtual lung epithelium within PhysiCell, an open source physics-based cell simulator ^[72]. In the project, epithelial cells are modelled as susceptible to SARS-CoV-2 infection. Other cell types, such as immune cells, are also included in the model but are not susceptible to the virus.

The PhysiCell SARS-CoV-2 framework consists of a hybrid agent-based model that aims to capture; virus transmission in the tissue, intracellular viral processes (e.g., binding, endocytosis, replication, exocytosis), infected cell responses (e.g., metabolism, chemokine secretion, cell death via apoptosis, cell death via pyroptosis), immune cell dynamics (e.g., recruitment, activation, infiltration and predation of infected cells) and tissue damage due to infection and/or host responses. For more details on each model component, we refer the reader to the most recent preprint from the coalition ^[55]. In line with a rapid prototyping and open source research approach, the project has been through several versions and will be updated further as each submodel is refined. Here, we will give a brief overview of the hybrid framework and discuss some of the implications of pyroptosis in the PhysiCell tissue simulator, using the version 3 code, with our added pyroptosis model. An initial iteration of our pyroptosis model is integrated into the version 4 tissue simulator code, and the refined model considered here will be implemented into future versions of the code.

3.2.   Simulation outline and implementation

At the start of each PhysiCell SARS-CoV-2 tissue simulation, virions (virus particles) are non-uniformly distributed across the epithelial cells. Virions diffuse across the tissue and can be endocytosed via ACE-2 receptors, and are thereafter replicated and released by epithelial cells. If the viral RNA in an epithelial cell surpasses a threshold value $R_{\rm{RNA}}$, the epithelial cell undergoes pyroptosis with a probability $\Pi \in [0, 1]$, and apoptosis with a probability $1-\Pi$. The intracellular pyroptosis model governs the cellular pyroptosis-induced secretion of the cytokines IL-18 and IL-1$\beta$, as well as the increase in cytoplasmic volume due to water influx, and the ultimate cell membrane rupture.

Once secreted, extracellular IL-18 and IL-1$\beta$ are modelled as diffusible fields across the tissue. IL-1$\beta$ has been shown to causes a pyroptotic bystander effect, whereby increased IL-1$\beta$ levels can induce pyroptosis in uninfected, or infected, cells ^{[13,17,73,74]}. Therefore, in the model, epithelial cells can internalise extracellular IL-1$\beta$, and if the level of internalised IL-1$\beta$ surpasses a threshold value $R_{\rm{IL-1}\beta}$, the pyroptosis cascade is initiated in the cell. The diffusion, secretion and uptake of virions and cytokines are handled by built-in PhysiCell functionalities. Furthermore, the PhysiCell framework has an inbuilt apoptosis mechanism which results in cell volume shrinkage until the cell is no longer viable. We include further details of this apoptosis model in Appendix E. The system of ODEs describing the pyroptosis pathway is numerically solved using the backward Euler method. Further details regarding implementation, as well as information on how to access the code, are provided in Appendix E.

3.3.   Agent-based model results

We now utilise the PhysiCell SARS-CoV-2 tissue simulator to investigate the effects of pyroptosis on a cellular-level scale. We design in silico experiments to study research questions related to virus-induced cell death. In these in silico experiments, we investigate scenarios in which epithelial cells are exposed to virus in the absence of an immune system, thus simulating in vitro experiments. We investigate the impact of two key parameters, $R_{\rm{IL-1}\beta}$ and $\Pi$, which are described in the above section.

Simulation results are provided in form of cell maps and cell count data in Figures 8 and 9. In the cell maps provided throughout this section, the results of a single run of the simulation are shown. Epithelial cells are shown in blue when they are viable, i.e., neither pyroptosing nor apoptosing. Uninfected epithelial cells are shown in blue, whilst infected epithelial cells are depicted as darker blue proportional to their viral load. Cells in which the pyroptosis process has been initiated due to viral loads or bystander effects (initiated by internalised IL-1$\beta$ from the external environment) are respectively shown in orange and red. Once cell membrane rupture occurs, as determined by the intracellular pyroptosis pathway model, a pyroptosing cell is removed from the system. Apoptosing cells turn black and reduce in size until they are removed from the system. The simulation data shown in the time-course plots is the total number of cells that are alive, and the accumulative number of cells that have initiated the investigated forms of cell death. In these time-course plots, the solid coloured lines denote the mean value of 10 runs of the simulation, and the standard deviation is plotted as a more transparent shade of the same colour. For all cases, the standard deviation between runs is generally small and thus, for clarity, we include a 'zoomed-in' version of each of the time-course plots in Appendix E, as seen in Figure E.2.

We first study the impact of the bystander effect of IL-1$\beta$, where pyroptosis may be induced in infected or uninfected cells that are close to pyroptosing cells that are secreting the cytokine. In Figure 8, we display the results corresponding to three values of the threshold, $R_{\rm{IL-1}\beta}$, specifically (A) $R_{\rm{IL-1}\beta} = 0.5\ R_{\rm{est}}$, (B) $R_{\rm{IL-1}\beta} = R_{\rm{est}}$ and (C) $R_{\rm{IL-1}\beta} = 1.5\ R_{\rm{est}}$, where $R_{\rm{est}}$ denotes an chosen threshold concentration. Comparing the three cell map progressions at the 30 h time-point, highlights qualitative differences between the results of the three cases. At this time-point, when we have a lower value of $R_{\rm{IL-1}\beta}$, as seen in case (A), we observe a larger number of pyroptosing cells (mostly induced by the bystander effect) and fewer viable cells in the system, in comparison to the other cases, (B) and (C). In the median case (B), there is some bystander effect occurring, however not to as large an extent as case (A). When we have a higher value of $R_{\rm{IL-1}\beta}$, as seen in case (C), we observe a lower bystander effect at the 30 h time-point. These qualitative differences are continued when comparing the 36 h time-point, whereby a lower value of $R_{\rm{IL-1}\beta}$ (case (A)), results in all epithelial cells being removed from the system, while a higher value of $R_{\rm{IL-1}\beta}$ (case (C)), results in a larger number of both viable and pyroptosing epithelial cells remaining. When comparing the time-course plots, we see that there is a significant difference in the number of cells that are induced to pyroptose through the bystander effect over time, where the bystander effect is the largest in case (A), and the smallest in case (C), as expected. We also observe that the time it takes for all viable cells to be removed from the system increases as we increase the threshold concentration required to induce pyroptosis through the bystander effect. These results highlight that the IL-1$\beta$ threshold concentration $R_{\rm{IL-1}\beta}$ affects the ratio of the method of pyroptosis inducement (virus or IL-1$\beta$), and also the time it takes for all cells to be removed from the system, as pyroptosis is a more rapid mode of cell death than apoptosis. Although the results from these simulations are qualitative, the model provides a framework in which to study bystander effects. Upon availability, data related to the IL-1$\beta$ levels that are required for pyroptosis induction could be incorporated in the model for a more quantitative analysis.

We next study the impact of the probability that an infected cell undergoes pyroptosis, as opposed to apoptosis, in response to viral load, as determined by the parameter $\Pi$. In Figure 9 we display the results corresponding to three values of $\Pi$. As before, comparing the three cell map progressions, the qualitative differences between the three cases are clear when comparing the cell maps at the 36 h time-point. In case (A) where the probability of pyroptosis is $\Pi = 0$, the results show that a large number of viable epithelial cells remain in the system at 36 h. We observe that there are fewer viable cells and more pyroptosing cells in case (B) than in case (A) at 36 h, moreover the bystander effect is the main cause of pyroptosis in case (B). Finally, when $\Pi = 1$ in case (C), by the 36 h time-point, all of the viable epithelial cells have been removed from the system, and only pyrotosing cells remain. Through the time-course plots, we can see a clear distinction between the number of cells undergoing apoptosis versus pyroptosis in the three cases. These simulation results show that increasing the value of $\Pi$ speeds up the overall eradication of epithelial cells. The reason for this is two-fold: pyroptosis contributes to bystander induced cell deaths, and the pyroptosis process is faster than apoptosis. More specifically, the time it takes for a cell to pyroptose or apoptose in the simulation is approximately 2 and 8.5 hours respectively. These results also highlight the impact of the bystander effect, which is, as expected, increased as $\Pi$ is increased, due to a larger number of viral load induced pyroptosing cells secreting IL-1$\beta$ into the external environment. As more cells are induced to pyroptose through the bystander effect, this further amplifies the levels of IL-1$\beta$ in the external environment, leading to an increased bystander effect. We additionally note, that in the time-course plot for case (B), where $\Pi = 0.5$, we observe a higher standard deviation resulting from 10 simulation runs. This is due to the added stochasticity in the model arising from the non-deterministic mode of cell death when $\Pi \neq 0$ and $\Pi \neq 1$. The time-course plots for all cases are included in a 'zoomed-in' format in Appendix E, Figure E.2.

4.   Conclusions and future work

Pyroptosis has been identified as a key mechanism involved in the cytokine storm and the inflammation associated with severe cases of COVID-19. Consequently, pyroptosis has been suggested as a target pathway to treat symptoms of COVID-19. In order to investigate the pathway of pyroptosis in further detail, we formulate a single-cell mathematical model that captures the key proteins and intracellular processes involved in this pro-inflammatory mode of cell death. Specifically, we model the process from DAMP/PAMP-induced TLR signalling to membrane rupture and the resulting cell death. The model is described in terms of a system of ODEs that captures the dynamics of intracellular molecule concentrations and cell volume. We further expand the model to include pharmacodynamic effects of a generic, NLRP3-targeting anti-inflammatory drug, and then perform both local and global sensitivity analyses in order to investigate how sensitive the model results are to the model parameters at their estimated values. Finally, we provide some detail on how this single-cell model is integrated into a wider cellular-scale model that is being used to investigate SARS-CoV-2 progression in lung epithelial cells ^[55].

Recently, it has been found that if the pyroptosis pathway is inhibited, a cell can instead undergo apoptosis, a form of non-inflammatory cell death ^{[11,22,24,26,75,76]}. In our single-cell model, described in Section 2, we have shown in our simulations that the inclusion of an NLRP3-targeting anti-inflammatory drug can delay the formation of the inflammasome, and therefore may provide a method of switching the cell death mode from inflammatory to non-inflammatory. The single-cell model developed in this study could further be extended to investigate more specific pharmacodynamic actions of anti-inflammatory drugs that target the formation of the NLRP3 inflammasome. A comprehensive review of such drugs is provided by Bertinara et al. ^[21]. One such drug is tranilast, an NLRP3-targeting drug that is traditionally used for treating inflammatory diseases, such as asthma, in Japan and South Korea. Tranislast is currently being re-purposed for treating COVID-19 related symptoms in a clinical trial^{* [2,77]}. Moreover, corticosteroids such as dexamethasone have been suggested as agents to counter the cytokine storm in severe cases of COVID-19 ^[78], and preliminary reports from clinical trials suggest that low-dose dexamethasone therapy has beneficial effects in COVID-19 patients requiring ventilation or supplemental oxygen ^[79]. As a corticosteroid, dexamethasone has many potential pharmacodynamic effects. Of relevance to the study of pyroptosis is that dexamethasone has been shown to directly bind to, and/or indirectly inhibit, NF-$\kappa$B, thus preventing the transcription of NLRP3 and pro-inflammatory cytokines ^[80]. Therefore, in future work, we could mathematically investigate the impact that this type of drug action would have on the single-cell pathway model, and we could further investigate the interplay between apoptosis and pyroptosis.

^*Chinese Clinical Trial Registry, registration number: ChiCDrug2000030002.

When investigating the role of pyroptosis in the PhysiCell SARS-CoV-2 framework, as described in Section 3, our results highlight that the mode of cell death impacts the cellular-scale dynamics. In this paper we have investigated an in vitro situation, however the PhysiCell framework supports the inclusion of within-host components, such as immune cells, to simulate in vivo scenarios. Previous studies have shown that IL-18 can act as a recruiter of immune cells ^{[13,14,15,16,17,18]}. Therefore, processes that allow extracellular IL-18 to act as an immune-cell recruiting cytokine can be implemented in the PhysiCell SARS-CoV-2 framework, when simulating in vivo scenarios.

Our model currently does not consider the effects of negative-feedback within the pyroptosis pathway. However, we note that there are several recent experimental studies that highlight the potential existence of negative-feedback loops that regulate the pyroptosis pathway at a transcriptional level ^[81,82]. For example, interferon regulatory factor 2 (IRF2) can bind to the pro-caspase-1 promoter, preventing pro-caspase-1 transcription, while IRF1 can upregulate this process ^[81]. Additionally, melatonin has been shown to inhibit NLRP3 priming ^[83,84]. Our model could be adapted to include some of these regulatory mechanisms, to further investigate the single-cell level dynamics of pyroptosis.

Investigating the role of interferons and interleukins within both the intracellular and cellular level would also be a beneficial extension to the model. At the intracellular level, it has been shown that interferon-1 (IFN-1) can decrease IL-18 transcription, and that interleukin-10 (IL-10) negatively regulates IL-1$\beta$ transcription ^[81]. Thus, at the cellular level, IFNs can play a role in inhibiting pyroptosis initiation, reducing the bystander effect and preventing cytokine storms ^[2,85,86]. The role of type Ⅰ IFNs in the treatment of COVID-19 is being investigated ^[87], therefore it would be beneficial to investigate these mechanisms in both the pathway model and in the wider PhysiCell multiscale model, where other mechanisms of IFNs are currently being considered.

We finally remark that pyroptosis has been identified as a potential contributor to symptoms in various diseases other than COVID-19 ^[2], such as other virus-induced diseases ^[17], auto-immune diseases ^[13,21] and cancer ^[23,26,76]. Therefore, the mathematical pyroptosis model developed in this paper could be modified and used to investigate the effect of pyroptosis in other disease scenarios.

Acknowledgments

SJH is supported by the Medical Research Council grant MR/R017506/1.

The authors would like to thank the SARS-CoV-2 Tissue Simulation Coalition for ongoing collaboration and feedback on model development. Further, this work is part of the RAMP (Rapid Assistance in Modelling the Pandemic) 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/.

Conflict of interest

The authors declare that there is no conflict of interest.

A.   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

● IFN: interferon

● IL-1$\beta$: interleukin 1$\beta$

● IL-18: interleukin 18

● LLR: leucine-rich repeat

● NBD: nucleotide-binding and oligomerisation domain

● NF-$\kappa$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$\beta$: pro-form of interleukin 1$\beta$

● pro-IL-18: pro-form of interleukin 18

● PYD: pyrin domain

● SARS-CoV-2: severe acute respiratory syndrome coronavirus 2

● TLR: toll like receptor

B.   Detailed model description and motivation

In this appendix we provide further details of the modelling terms chosen within our single-cell pathway model described in Section 2. Model variables/components are listed in Table B.1 and the model parameters are discussed in Appendix C.

B.1.   Signal 1.

The initial signal, $S_{1}$, as defined in Eq (2.2), is here considered to be either on, $S_1 = 1$, or off, $S_1 = 0$. This binary signal is motivated by the modelling assumption that there will be no inflammasome-inducing 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 assume that once $S_1$ is turned on, this signal is kept at a constant level until the cell dies.

B.2.   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. ^[88] showed that increasing the time that a cell was exposed to ATP (yielding $S_2$ to turn on) did not alter the subsequent inflammasome dynamics, nor the expression of involved proteins. Thus, assuming that $S_2$ is binary, i.e., on or off, throughout the simulation is a reasonable model assumption.

B.3.   Reaction terms.

To describe the transcription or cleavage terms we use a Hill function ^[59] of the form,

Here, A₅₀ represents the level of [A] required for the Hill function to reach the value of a half, often referred to the half-max value of [A] ^[89]. The rate parameter $\alpha$ determines the cleavage or transcription rate. To describe the activation terms in the ODE system we consider chemical reactions of the form,

where $A$ and $B$ are the inactive and active components, respectively. The reaction rate $k$ determines how fast this process occurs. Using the law of mass action we can then translate this into the term used in the ODE for concentration $[A]$, i.e., $- k[A]$, and for concentration $[B]$, i.e., $+ k[A]$.

Similarly, to include decay terms, we consider a chemical reaction of the form,

where $\emptyset$ here denotes that the component is removed from the system. Again, using the law of mass action we can then use this to write the decay terms within the ODE for concentration $[A]$, i.e., $- \delta [A]$.

Catalytic interactions within the model are added considering a reaction of the form,

Here, component $A$ and $B$ bind together briefly to transform $A$ into a 3rd component $C$. Using the law of mass action, we can then use this to write the reactions within the ODE for concentration $[A]$, i.e., $- k [A][B]$, and concentration $[C]$, i.e., $+k [A][B]$. Note, as $B$ is not used up the terms for concentration $[B]$ will will not be included.

We consider the NLRP3 oligomerisation process to be accumulative, that is once two active NLRP3 proteins bind, then any additional proteins will bind on to this bound NLRP3 forming the inflammasome. If p proteins are required to form the inflammasome the interactions will be of the form,

where $m$ is the number proteins bound together. Note, we have that [NLRP3₁] = [NLRP3_a] as this is unbound active protein.

The law of mass action gives for m = 2, $\ldots$, p-1,

and

The total level of bound or oligomerised protein is,

Considering the rate of change of [NLRP3_o] we obtain,

This simplifies to,

as all other terms cancel out.

B.4.   Anti-inflammatory drug terms.

The influence of the anti-inflammatory drug is included in the modified version of the model which includes additional terms for the dynamics of NLRP3_a, as given in Eq (2.27). Also included in this modified model are the dynamics of the anti-inflammatory drug, Drug, and the complex Drug $\cdot$ NLRP3_a, as listed in Eqs (2.28) and (2.29). Here, we consider a binding reaction of the form,

Mass action kinetics allows us to write this interaction in terms of ODEs. We consider the anti-inflammatory drug to be a covalent inhibitor, and thus we could add an extra step in the above reaction. However, in lieu of relevant data, we are currently considering only one reaction step for Drug binding to NLRP3_o, in an effort to minimise the model.

C.   Model parameter values

In this appendix, we provide further details on how the model parameter values, listed in Table C.1, were estimated.

C.1.   Decay parameters.

In ^[61], it was found that NLRP3 had a half-life of approximately 6 hrs when cells are stimulated with lipopolysaccharide (LPS), which translates to a decay rate of approximately 0.002 $\mu\rm{m} \rm{ min}^{-1}$. Therefore we set,

as we assume measure concentrations as arbitrary units. In ^[67], it was found that pro-IL-1$\beta$ had a half-life of approximately 2.5 hrs in human monocytes, which translates to a decay rate of approximately 0.004 $\mu\rm{m}\ \rm{min}^{-1}$. Therefore we set,

as we assume measure concentrations as arbitrary units.

C.2.   The timeline of pyroptosis

We use data from three experimental studies ^[56,60,68] to estimate a timeline for events regarded in our mathematical model. The values of the model parameters $\alpha_{1}-\alpha_{5}, \ \rm{C1}_{50}, \ \gamma_{\rm{C1}}, \ \gamma_{\rm{NF}}, \ k_{1}-k_{7}, \ \rm{NF}_{50}$ and $V_{c}$ are chosen to obtain model dynamics that follow this experimentally motivated timeline of events A-D, which is illustrated in Figure C.1. Further detail of constructing this timeline is provided in the following paragraphs.

Bagaev et al. ^[56] studied NF-$\kappa$B kinetics in primary macrophages subjected to bacterial lipopolysaccharide (LPS). de Vasconcelos et al. ^[60] investigated subcellular key events that define pyroptosis using single-cell analysis. 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. Martín-Sánchez et al. ^[68] measured the levels of pro-IL-1$\beta$ and external IL-1$\beta$ over time after inflammasome activation in systems with mouse BMDMs.

A: Bagaev et al. ^[56] reported that the nuclear NF-$\kappa$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.

B: de Vasconcelos et al. ^[60] showed that the influx of Ca$^{2+}$ occurs before total membrane permeabilisation in pyroptosis (see Figure 7 in de Vasconcelos et al. ^[60]). Further data from their study 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, as well as Supplementary Figure 3 in de Vasconcelos et al. ^[60]). 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 the time before cell rupture (end time) when the inflammasome forms (i.e. when the concentration of NLRP3_o reaches the threshold value $n$). Therefore event B occurs around end time - 43 minutes = (120-43) minutes = 77 minutes.

C: In de Vasconcelos et al.'s study, cell volumes increase gradually for approximately 13 minutes prior to membrane rupture. Furthermore, cell volumes were reported to increase by approximately 50%. (See Figure 6a, b in de Vasconcelos et al. ^[60]).

D: Supernatants from cell cultures were assayed for LDH activity in de Vasconcelos et al.'s study. 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. ^[60]). From this we approximate the time between signal 1 turning on in the model, and complete membrane rupture, to be approximately 120 minutes.

Furthermore, results from Martín-Sánchez et al.'s ^[68] study highlight that the release of IL-1$\beta$ coincided with membrane permeability (the formation of GSDMD-N derived membrane pores), and that eventually all of the pro-IL-1$\beta$ present at the start of the experiments was cleaved and released from the cell, with approximately 90% released within 120 minutes.

C.3.   Anti-inflammatory drug parameters.

We consider drug concentrations to be constant during the simulated time span, effectively describing an in vitro scenario without drug wash out. This could be altered to include drug dosages varying in time if appropriate. Tranilast and oridonin examples of anti-inflammatory drugs that can inhibit the formation of the NLRP3 inflammasome by binding to active NLRP3 ^[21]. Immunoblot data presented by Huang et al. ^[58], pertaining to LPS-primed bone marrow derived macrophages (BMDMs) treated with tranilast, demonstrate dose-dependent (25-100$\mu$M) drug responses, where the maximal inhibition was reached at at 100 $\mu$M of the drug. In our model, we implicitly assume that pyroptosis is inhibited if the membrane rupture does not occur within 300 simulated minutes. We here set $k_{-\rm{D}}/k_{+\rm{D}} = 0.01 \;{\rm{ a.u.}}$, so that pyroptosis is inhibited for drug concentrations over 0.75 a.u..

D.   MATLAB code and set up of numerical simulations for the single-cell model

The single-cell mathematical model was implemented in MATLAB using the ODE solver ode15s. The simulation results provided in Figures 3 and 4 are obtained by running the code with the initial conditions given in Section 2.9, along with the parameters given in Table C.1 and the signals $S_1$ and $S_2$ turned on. The MATLAB code used to run these simulations is available on our project GitHub page: https://github.com/Pyroptosis.

E.   Details of the PhysiCell implementation

To implement the ODE pathway model, as described in Section 2, into the PhysiCell framework, we discretise the full system (2.31) with functions (2.32) using the Backward Euler method. Note, that we do not consider the drug in the PhysiCell setting, that is we let $[\rm{Drug}](t)\equiv 0$. The parameter values used are those previously shown in Table C.1. The base threshold level of IL-1$\beta$ required to initiate pyroptosis is denoted by parameter $R_{\rm{est}}$ and the threshold level of viral RNA required to initiate viral death is denoted by $R_{\rm{RNA}}$. For the results shown in Figures 8 and 9, we set,

The values of the parameters $\Pi$ and $R_{\rm{IL-1}\beta}$ are varied, and shown in the caption of each Figure. We additionally plot the solutions of the Backward Euler method using these parameters in Figure E.1. For more details of the full PhysiCell model, we refer the reader to the most recent preprint ^[55], or the SARS-CoV-2 Tissue Simulation Coalition webpage: http://physicell.org/covid19/. Additionally, the MATLAB code used to create Figure E.1, and the specific PhysiCell C++ code used to create Figures 8 and 9 are available on our project GitHub page: https://github.com/Pyroptosis.

E.1.   Apoptosis model in PhysiCell

The following details are taken from ^[72], and describes the in built standard apoptosis model in PhysiCell. Cell volume is tracked through multiple properties including,

● $V_F$ - total fluid volume of the cell,

● $V_{CS}$ - total solid cytoplasmic volume of the cell,

● $V_{NS}$ - total solid nuclear volume of the cell.

Once apoptosis initiates, the cell volume changes according to the following equations,

where the parameters $r_F, \ r_{CS}, \ r_{NS}$ are the rate of change of the total fluid volume, the solid cytoplasmic volume and the solid nuclear volume, respectively.

In the PhysiCell model, the cell is removed once a threshold volume $V_m$ has been reached, or the cell has been apoptosing for longer than time $T_A$. The default parameter values for apoptosis of epithelial cells ^[72] are,

We consider that along with natural cell death, that the virus can induce apoptosis or pyroptosis once viral RNA has reached a certain threshold, $R_{\rm{RNA}}$.