A single-cell mathematical model of SARS-CoV-2 induced pyroptosis and the
effects of anti-inflammatory intervention

Sara J Hamis,1 Fiona R Macfarlane.1

1 School of Mathematics and Statistics, University of St Andrews, St Andrews,
Scotland, UK.

Both authors contributed equally to this work.

This is version 4 of the pre-print article.

###### Abstract

Pyroptosis is an inflammatory mode of cell death that can contribute to the
cytokine storm associated with severe cases of coronavirus disease 2019
(COVID-19). The formation of the NLRP3 inflammasome is central to pyroptosis,
which may be induced by severe acute respiratory syndrome coronavirus 2 (SARS-
CoV-2). Inflammasome formation, and by extension pyroptosis, may be inhibited
by certain anti-inflammatory drugs.

In this study, we present a single-cell mathematical model that captures the
formation of the NLRP3 inflammasome, pyroptotic cell death and responses to
anti-inflammatory intervention that hinder the formation of the NLRP3
inflammasome. The model is formulated in terms of a system of ordinary
differential equations (ODEs) that describe the dynamics of the key biological
components involved in pyroptosis. Our results demonstrate that an anti-
inflammatory drug can delay the formation of the NLRP3 inflammasome, and thus
may alter the mode of cell death from inflammatory (pyroptosis) to non-
inflammatory (e.g., apoptosis). The single-cell model is implemented within a
SARS-CoV-2 tissue simulator, in collaboration with a multidisciplinary
coalition investigating within host-dynamics of COVID-19. In this paper, we
additionally provide an overview of the SARS-CoV-2 tissue simulator and
highlight the effects of pyroptosis on a cellular level.

Keywords: Pyroptosis; COVID-19; NLRP3 Inflammasome; Cytokine Storm;
Mathematical Modelling.

##  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, 19, 20, 9], 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 [4, 2, 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 [24, 17, 25, 21, 26].

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 model components. Blue arrows represent
facilitation of component formation or reactions. Red arrows represent signals
that can be turned on or off in the model. GSDMD-N induced membrane pores
(shown in red) allow for the outflux of inflammatory cytokines and the influx
of extracellular water, causing the cell to swell and the cell membrane to
eventually rupture. Our mathematical model includes the effect of a drug that
inhibits the NLRP3o inflammasome base from forming. Figure abbreviations are
listed in Appendix A.

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

Figure 2: NLRP3 inflammasome formation. (a) Through a wide array of stimuli,
inactive NLRP3 can become activated. Active NLRP3 can bind together to form
the inflammasome base, this binding can be inhibited by NLRP3-targeting anti-
inflammatory drugs. (b) Once the inflammasome base is formed, the pyrin domain
(PYD) of ASC molecules can bind to the PYD of NLRP3 molecules in the
inflammasome base. Thereafter, the caspase activation and recruitment domain
(CARD) of pro-caspase-1 can bind to the CARD of NLRP3-bound ASC molecules.
This results in the dimerisation and activation of caspase-1. Abbreviations
used in the figure are listed in Appendix 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,

$\text{[NF-$\kappa$B\textsubscript{n}]}({t})=\text{[NF-$\kappa$B\textsubscript{n}]}(0)+S_{1}\cdot
h\ e^{-\frac{\log^{2}(t/\tau)}{s}},$ (1)

where [NF-$\kappa$Bn](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,

$S_{1}:=\begin{cases}1&\mbox{{\it on}, if a DAMP/PAMP induced signal is
present, }\\\ 0&\mbox{{\it off,} otherwise.}\end{cases}$ (2)

In Equation (1), $h$ is a constant that denotes the peak elevation of the
[NF-$\kappa$Bn] 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$Bn] peak over time. Accordingly, the parameters $h$, $\tau$ and
$s$ can be varied to adjust the shape of the [NF-$\kappa$Bn] peak, as is shown
in the Supplementary Material (Supplementary Material, S1).

###  2.2 NLRP3 dynamics

Following nuclear NF-$\kappa$B translocation, inactive NLRP3, NLRP3i, 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_{\text{NF}}$ and half-max concentration-value NF50. 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, NLRP3a molecules can
oligomerise (bind together) to form a wheel-shaped inflammasome base, modelled
by the concentration of oligomerised NLRP3, i.e., NLRP3o. In the model, we
thus consider NLRP3 protein concentrations in three different forms: [NLRP3i],
[NLRP3a] and [NLRP3o]. Inactive and active NLRP3 decays in the model at a rate
$\delta_{1}$, however we do not include the decay of NLRP3o in the model. The
reason behind this choice of omitting NLRP3o decay is that, in experimental
works, the functions of the NLRP3 inflammasome are evident until cell death in
the context of pyroptosis [62, 61, 60]. Furthermore, it is unclear whether
disassembly of the inflammasome is possible in any context [64, 63]. 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
[NLRP3i], [NLRP3a] and [NLRP3o] over time as,

$\displaystyle\frac{d\text{[NLRP3\textsubscript{i}]}}{dt}$ $\displaystyle=$
$\displaystyle\alpha_{1}\frac{[\overline{\text{NF-$\kappa$B\textsubscript{n}}}]^{\gamma_{\text{NF}}}}{\text{NF}_{50}^{\gamma_{\text{NF}}}+[\overline{\text{NF-$\kappa$B\textsubscript{n}}}]^{\gamma_{\text{NF}}}}-S_{2}\
k_{1}\ \text{[NLRP3\textsubscript{i}]}-\delta_{1}\
\text{[NLRP3\textsubscript{i}]},$ (3)
$\displaystyle\frac{d\text{[NLRP3\textsubscript{a}]}}{dt}$ $\displaystyle=$
$\displaystyle S_{2}\ k_{1}\ \text{[NLRP3\textsubscript{i}]}-k_{2}\
\text{[NLRP3\textsubscript{a}]}^{2}-\delta_{1}\
\text{[NLRP3\textsubscript{a}]},$ (4)
$\displaystyle\frac{d\text{[NLRP3\textsubscript{o}]}}{dt}$ $\displaystyle=$
$\displaystyle k_{2}\ \text{[NLRP3\textsubscript{a}]}^{2},$ (5)

where [$\overline{\text{NF-$\kappa$B\textsubscript{n}}}$] denotes the
deviation from baseline in [NF-$\kappa$Bn] so that,

$[\overline{\text{NF-$\kappa$B\textsubscript{n}}}](t)=\text{[NF-$\kappa$B\textsubscript{n}]}(t)-\text{[NF-$\kappa$B\textsubscript{n}]}(0).$
(6)

The binary signal $S_{2}$ is set as,

$S_{2}:=\begin{cases}1&\mbox{{\it on}, if an activation signal is present,
}\\\ 0&\mbox{{\it off,} otherwise.}\end{cases}$ (7)

###  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 [ASCf] and [ASCb], respectively. We
consider this binding process to be irreversible and to occur at a rate
$k_{3}$, hence the rate of change of [ASCf] and [ASCb] are here described by,

$\displaystyle\frac{d\text{[ASC\textsubscript{f}]}}{dt}$ $\displaystyle=$
$\displaystyle-k_{3}\ F([\text{NLRP3\textsubscript{o}}])\
\text{[NLRP3\textsubscript{o}]}\ \text{[ASC\textsubscript{f}]},$ (8)
$\displaystyle\frac{d\text{[ASC\textsubscript{b}]}}{dt}$ $\displaystyle=$
$\displaystyle k_{3}\ F([\text{NLRP3\textsubscript{o}}])\
\text{[NLRP3\textsubscript{o}]}\ \text{[ASC\textsubscript{f}]},$ (9)

where the function $F([\text{NLRP3\textsubscript{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,

$F(\text{[NLRP3\textsubscript{o}])}=\frac{1}{1+\Big{(}\frac{\text{[NLRP3\textsubscript{o}]}+a}{b}\Big{)}^{-c}}.$
(10)

Here, values $a=1$ (a.u.) and $b=2$ (a.u.) can be used to obtain
$F([\text{NLRP3\textsubscript{o}}])\approx 0$ for
$[\text{NLRP3\textsubscript{o}}]<1$ and
$F(\text{[}\text{NLRP3\textsubscript{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
$[\text{NLRP3\textsubscript{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,

$\text{[ASC\textsubscript{f}]}+\text{[ASC\textsubscript{b}]}=\text{constant}.$
(11)

Free ASC levels are adequate for pyroptosis at homeostasis [24, 21].
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,

$\displaystyle\frac{d\text{[pro-C1]}}{dt}$ $\displaystyle=$ $\displaystyle-
k_{4}\ \text{[ASC\textsubscript{b}]}\ \text{[pro-C1]},$ (12)
$\displaystyle\frac{d\text{[C1]}}{dt}$ $\displaystyle=$ $\displaystyle k_{4}\
\text{[ASC\textsubscript{b}]}\ \text{[pro-C1]},$ (13)

so that,

$\text{[pro-C1]}+\text{[C1]}=\text{constant}.$ (14)

Pro-caspase-1 is present in adequate levels in the cell prior to cellular
responses to DAMPs or PAMPs [24, 21]. 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_{\text{C1}}$ is the Hill coefficient, and C150 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,

$\displaystyle\frac{d\text{[GSDMD]}}{dt}$ $\displaystyle=$
$\displaystyle-\alpha_{2}\
\frac{\text{[C1]}^{\gamma_{C1}}}{\text{C1}_{50}^{\gamma_{C1}}+\text{[C1]}^{\gamma_{C1}}}\
\text{[GSDMD]},$ (15) $\displaystyle\frac{d\text{[GSDMD-N]}}{dt}$
$\displaystyle=$ $\displaystyle\alpha_{2}\
\frac{\text{[C1]}^{\gamma_{C1}}}{\text{C1}_{50}^{\gamma_{C1}}+\text{[C1]}^{\gamma_{C1}}}\
\text{[GSDMD]},$ (16)

so that the following conservation law holds,

$\text{[GSDMD]}+\text{[GSDMD-N]}=\text{constant}.$ (17)

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/18c] and [IL-1$\beta$e/18e]. 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-18c 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$Bn], in a similar way to NLRP3, therefore we describe this process
using a Hill function with constant coefficient $\alpha_{3}$, Hill coefficient
$\gamma_{\text{NF}}$ and the half-max concentration of NF-$\kappa$B, NF50. 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,

$\displaystyle\frac{d\text{[pro-IL-1}\beta]}{dt}$ $\displaystyle=$
$\displaystyle\alpha_{3}\
\frac{[\overline{\text{NF-$\kappa$B\textsubscript{n}}}]^{\gamma_{NF}}}{\text{NF}_{50}^{\gamma_{NF}}+[\overline{\text{NF-$\kappa$B\textsubscript{n}}}]^{\gamma_{NF}}}-\alpha_{4}\
\frac{\text{[C1]}^{\gamma_{C1}}}{\text{C1}_{50}^{\gamma_{C1}}+\text{[C1]}^{\gamma_{C1}}}\
\text{[pro-IL-1$\beta$]}$ (18) $\displaystyle-\delta_{2}\ \text{[pro-
IL-1$\beta$]},$
$\displaystyle\frac{d\text{[IL-1$\beta$\textsubscript{c}]}}{dt}$
$\displaystyle=$ $\displaystyle\alpha_{4}\
\frac{\text{[C1]}^{\gamma_{C1}}}{\text{C1}_{50}^{\gamma_{C1}}+\text{[C1]}^{\gamma_{C1}}}\
\text{[pro-IL-1$\beta$]}-k_{5}\ G\ \text{[IL-1$\beta$\textsubscript{c}]}$ (19)
$\displaystyle-\delta_{2}\ \text{[IL-1$\beta$\textsubscript{c}]},$
$\displaystyle\frac{d\text{[IL-1$\beta$\textsubscript{e}]}}{dt}$
$\displaystyle=$ $\displaystyle k_{5}\ G\text{
[IL-1$\beta$\textsubscript{c}]},$ (20)

where $[\overline{\text{NF-$\kappa$B\textsubscript{n}}}]$ is defined in
Equation (6) and $G$ is a non-dimensionalised value that allows for
[GSDMD-N]-dependent cytokine outflux, and water influx, such that,

$G=G(\text{[GSDMD-N]}):=\frac{\text{[GSDMD-N]}}{\text{[GSDMD]+[GSDMD-N]}}.$
(21)

##### IL-18:

We describe the rate of change of [pro-IL-18], [IL-18c] and [IL-18e] as,

$\displaystyle\frac{d\text{[pro-IL-18]}}{dt}$ $\displaystyle=$
$\displaystyle-\alpha_{5}\
\frac{\text{[C1]}^{\gamma_{C1}}}{\text{C1}_{50}^{\gamma_{C1}}+\text{[C1]}^{\gamma_{C1}}}\
\text{[pro-IL-18]},$ (22)
$\displaystyle\frac{d\text{[IL-18\textsubscript{c}]}}{dt}$ $\displaystyle=$
$\displaystyle\alpha_{5}\
\frac{\text{[C1]}^{\gamma_{C1}}}{\text{C1}_{50}^{\gamma_{C1}}+\text{[C1]}^{\gamma_{C1}}}\
\text{[pro-IL-18]}-k_{6}\ G\ \text{[IL-18\textsubscript{c}]},$ (23)
$\displaystyle\frac{d\text{[IL-18\textsubscript{e}]}}{dt}$ $\displaystyle=$
$\displaystyle k_{6}\ G\ \text{[IL-18\textsubscript{c}]},$ (24)

so that,

$\text{[pro-
IL-18]}+\text{[IL-18\textsubscript{c}]}+\text{[IL-18\textsubscript{e}]}=\text{constant}.$
(25)

The levels of pro-IL-18 in the cell are kept at a homeostatic level [24, 21]
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,

$\displaystyle\frac{dV}{dt}$ $\displaystyle=$ $\displaystyle k_{7}\ G\ V.$
(26)

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 NLRP3a
oligomerisation. To incorporate this in our model, we consider the drug
binding to NLRP3a to form the complex Drug $\cdot$ NLRP3a. The forward and
reverse rate constants for this reaction are denoted by $k_{+\text{D}}$ and
$k_{-\text{D}}$, respectively. When including these drug mechanisms in our
model, Equation (4) describing the dynamics of NLRP3a must be modified to
include two additional terms,

$\frac{d\text{[NLRP3\textsubscript{a}]}}{dt}=\ldots-k_{+\text{D}}\
\text{[Drug]}\ \text{[NLRP3\textsubscript{a}]}+k_{-\text{D}}\ \text{[Drug
$\cdot$ NLRP3\textsubscript{a}}].$ (27)

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

$\displaystyle\frac{d\text{[Drug]}}{dt}$ $\displaystyle=$ $\displaystyle-
k_{+\text{D}}\ \text{[Drug]}\ \text{[NLRP3\textsubscript{a}]}+k_{-\text{D}}\
\text{[Drug $\cdot$ NLRP3\textsubscript{a}]},$ (28)
$\displaystyle\frac{d\text{[Drug $\cdot$ NLRP3\textsubscript{a}]}}{dt}$
$\displaystyle=$ $\displaystyle k_{+\text{D}}\ \text{[Drug]}\
\text{[NLRP3\textsubscript{a}]}-k_{-\text{D}}\ \text{[Drug $\cdot$
NLRP3\textsubscript{a}]},$ (29)

where the following conservation law holds,

[Drug]+[Drug $\cdot$ NLRP3a] = constant. (30)

###  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 (11), (14), (17), (25) and (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,

$\displaystyle\text{[NF-$\kappa$B\textsubscript{n}]}({t})$ $\displaystyle=$
$\displaystyle\text{[NF-$\kappa$B\textsubscript{n}]}(0)+S_{1}\cdot h\
e^{-\frac{\log^{2}(t/\tau)}{s}},$
$\displaystyle\frac{d\text{[NLRP3\textsubscript{i}]}}{dt}$ $\displaystyle=$
$\displaystyle\alpha_{1}\text{Hill}_{\text{NF}}-S_{2}k_{1}\text{[NLRP3\textsubscript{i}]}-\delta_{1}\text{[NLRP3\textsubscript{i}]},$
$\displaystyle\frac{d\text{[NLRP3\textsubscript{a}]}}{dt}$ $\displaystyle=$
$\displaystyle
S_{2}k_{1}\text{[NLRP3\textsubscript{i}]}-k_{2}\text{[NLRP3\textsubscript{a}]}^{2}-\delta_{1}\text{[NLRP3\textsubscript{a}]}$
$\displaystyle-k_{+\text{D}}\ \text{[Drug]}\
\text{[NLRP3\textsubscript{a}]}+k_{-\text{D}}\ \text{[Drug $\cdot$
NLRP3\textsubscript{a}}],$
$\displaystyle\frac{d\text{[NLRP3\textsubscript{o}]}}{dt}$ $\displaystyle=$
$\displaystyle k_{2}\text{[NLRP3\textsubscript{a}]}^{2},$
$\displaystyle\frac{d\text{[ASC\textsubscript{b}]}}{dt}$ $\displaystyle=$
$\displaystyle k_{3}\ F_{\text{N}_{\text{o}}}\
\text{[NLRP3\textsubscript{o}]}\ (1-\text{[ASC\textsubscript{b}]}),$
$\displaystyle\frac{d\text{[C1]}}{dt}$ $\displaystyle=$ $\displaystyle k_{4}\
\text{[ASC\textsubscript{b}]}\ (1-\text{[C1]}),$
$\displaystyle\frac{d\text{[GSDMD-N]}}{dt}$ $\displaystyle=$
$\displaystyle\alpha_{2}\ \text{Hill}_{\text{C1}}\ (1-\text{[GSDMD-N]}),$
$\displaystyle\frac{d\text{[pro-IL-1}\beta]}{dt}$ $\displaystyle=$
$\displaystyle\alpha_{3}\ \text{Hill}_{\text{NF}}-\alpha_{4}\
\text{Hill}_{\text{C1}}\text{[pro-IL-1$\beta$]}-\delta_{2}\ \text{[pro-
IL-1$\beta$]},$ (31)
$\displaystyle\frac{d\text{[IL-1$\beta$\textsubscript{c}]}}{dt}$
$\displaystyle=$ $\displaystyle\alpha_{4}\ \text{Hill}_{\text{C1}}\text{[pro-
IL-1$\beta$]}-k_{5}G\text{[IL-1$\beta$\textsubscript{c}]}-\delta_{2}\
\text{[IL-1$\beta$\textsubscript{c}]},$
$\displaystyle\frac{d\text{[IL-1$\beta$\textsubscript{e}]}}{dt}$
$\displaystyle=$ $\displaystyle k_{5}\
G\text{[IL-1$\beta$\textsubscript{c}]},$
$\displaystyle\frac{d\text{[IL-18\textsubscript{c}]}}{dt}$ $\displaystyle=$
$\displaystyle\alpha_{5}\text{Hill}_{\text{C1}}(1-\text{[IL-18\textsubscript{c}]}-\text{[IL-18\textsubscript{e}]})-k_{6}G\text{[IL-18\textsubscript{c}]},$
$\displaystyle\frac{d\text{[IL-18\textsubscript{e}]}}{dt}$ $\displaystyle=$
$\displaystyle k_{6}G\text{[IL-18\textsubscript{c}]},$
$\displaystyle\frac{d\text{V}}{dt}$ $\displaystyle=$ $\displaystyle
k_{7}G\text{V},$ $\displaystyle\frac{d\text{[Drug]}}{dt}$ $\displaystyle=$
$\displaystyle-k_{+\text{D}}\ \text{[Drug]}\
\text{[NLRP3\textsubscript{a}]}+k_{-\text{D}}\ \text{[Drug $\cdot$
NLRP3\textsubscript{a}]},$ $\displaystyle\frac{d\text{[Drug $\cdot$
NLRP3\textsubscript{a}]}}{dt}$ $\displaystyle=$ $\displaystyle k_{+\text{D}}\
\text{[Drug]}\ \text{[NLRP3\textsubscript{a}]}-k_{-\text{D}}\ \text{[Drug
$\cdot$ NLRP3\textsubscript{a}]},$

with the following definitions,

$\displaystyle\text{Hill}_{\text{NF}}$ $\displaystyle=$
$\displaystyle\frac{(\text{[NF-$\kappa$B\textsubscript{n}]}(t)-\text{[NF-$\kappa$B\textsubscript{n}]}(0))^{\gamma_{\text{NF}}}}{\text{NF}_{50}^{\gamma_{\text{NF}}}+(\text{[NF-$\kappa$B\textsubscript{n}]}(t)-\text{[NF-$\kappa$B\textsubscript{n}]}(0))^{\gamma_{\text{NF}}}},$
$\displaystyle F_{\text{N}_{\text{o}}}$ $\displaystyle=$
$\displaystyle\frac{1}{1+\left(\frac{\text{[NLRP3\textsubscript{o}]}-a}{b}\right)^{-c}},$
(32) $\displaystyle\text{Hill}_{\text{C1}}$ $\displaystyle=$
$\displaystyle\frac{\text{[C1]}^{\gamma_{C1}}}{\text{C1}_{50}^{\gamma_{C1}}+\text{[C1]}^{\gamma_{C1}}}.$

The reduced model is implemented in MATLAB [66], where the system of ODEs,
(31) with (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,

$[\text{NF-$\kappa$B}_{\text{n}}](0)=0.25\text{ a.u.},$

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

$[\text{ASC\textsubscript{f}}](0)=1\text{ a.u.},\ [\text{pro-C1}](0)=1\text{
a.u.},$ $\ \text{[GSDMD]}(0)=1\text{ a.u.},\ \text{[pro-IL-18]}(0)=1\text{
a.u.},\ V(0)=1\text{ a.u.},$

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 (31) with functions (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$Bn] levels have peaked, the transcription, and
by extension the synthesis, of NLRP3i and pro-IL-1$\beta$ are reduced as
NF-$\kappa$B leaves the nucleus. Synthesised NLRP3i activates to become
NLRP3a, which in turn oligomerises and binds together to form the inflammasome
base, here expressed in terms of the concentration of NLRP3o. These dynamics
are captured in the subplot displaying the time evolution of [NLRP3i],
[NLRP3a] and [NLRP3o]. 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 [ASCb] levels
increase when [NLRP3o] reaches the value $n$, where here $n=1$ a.u..
Furthermore, as [ASCb] 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-18c]
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-18c] respectively decrease in favour of the
extracellular concentrations [IL-1$\beta$e] and [IL-18e]. 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.

Figure 3: Numerical simulations results of the system of ODEs (31) with
functions (32) and [Drug](0)=0 a.u., that is, the model of pyroptosis in the
absence of the anti-inflammatory drug. We display the concentration of each
model component in arbitrary units (a.u.) over time, as well as cell volume
dynamics. Note that the time progression stops once the critical volume
$V_{c}$ is reached and the ultimate cell membrane rupture occurs.

We next consider numerical simulations of the model including the effects of
an anti-inflammatory drug, that is system of ODEs (31) with functions (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 [NLRP3o] 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.

Figure 4: Numerical simulations results of the system of ODEs (31) with
functions (32) with non-zero drug initial concentrations. We display the
concentration of each model component in arbitrary units (a.u.) over time, as
well as cell volume dynamics for several dosages of the drug. The dynamics of
the key model components are plotted, where the line colour corresponds to
applied drug dosage, as described in the legend. Note that the time
progression stops if the critical volume (i.e., $V_{c}=1.5$ a.u.) is reached,
as the cell membrane ruptures.

###  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 (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 $\pm 20\%$ 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.

Figure 5: Robustness analyses. The plots show numerical simulations results of
the ODE system (31), where the input parameter $\alpha_{1}$ (in a.u. min-1) is
perturbed, whilst other model parameters are fixed at their estimated values.
Model components are measured in terms of arbitrary concentration units (a.u.)
or volume. Note that the time progression stops once the critical volume
$V_{c}$=1.5 a.u. is reached and the ultimate cell membrane rupture occurs.

####  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 NLRP3o concentration at 77 minutes
(the time point where the inflammasome base is formed i.e., when [NLRP3o]
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.

Figure 6: Global sensitivity analysis is performed using Latin hypercube
sampling and analysis. The output, i.e., the NLRP3o concentration (in
arbitrary units of concentration) at 77 minutes, is plotted over three input
parameters. Each input parameter range spans $\pm 10\%$ of the estimated
parameter value. The Pearson correlation coefficient (R), measuring the linear
association between each input-output pair, is provided at the top right
corner of each corresponding subplot.

####  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 (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 NLRP3a levels increase at the start of the simulation (for
times under $\sim$ 100 minutes). NLRP3o levels also increase with
$\alpha_{1}$, where the magnitude of this increase plateaus when the
inflammasome base is formed. Thus the outputs [NLRP3a] and [NLRP3o] 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$.

Figure 7: A derivative-based sensitivity analysis shows that various output
results are sensitive to perturbations of the input parameter $\alpha_{1}$,
for the estimated parameter values of $\Theta$. Here $z^{i}_{j}$ denotes the
sensitivity of the output concentration computed in equation $i$, with respect
to the input parameter $\theta_{j}$, where $\xi$(compound concentration) maps
each compound to its governing equation index $i$.

##  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 subgroups 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_{\text{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, 73, 74, 17]. 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_{\text{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_{\text{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_{\text{IL-1}\beta}$, specifically (A) $R_{\text{IL-1}\beta}=0.5\
R_{\text{est}}$, (B) $R_{\text{IL-1}\beta}=R_{\text{est}}$ and (C)
$R_{\text{IL-1}\beta}=1.5\ R_{\text{est}}$, where $R_{\text{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_{\text{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_{\text{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_{\text{IL-1}\beta}$ (case (A)), results in all epithelial cells being
removed from the system, while a higher value of $R_{\text{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_{\text{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.

Figure 8: Cell map progressions and averaged time-course data for three
situations, investigating the bystander effect of IL-1$\beta$ in the PhysiCell
SARS-CoV-2 framework. We test three values for the threshold level of
IL-1$\beta$ required to initiate pyroptosis: (A) $R_{\text{IL-1}\beta}=0.5\
R_{\text{est}}$, (B) $R_{\text{IL-1}\beta}=R_{\text{est}}$ and (C)
$R_{\text{IL-1}\beta}=1.5\ R_{\text{est}}$, where $R_{\text{est}}$ is a chosen
base-line value. Viable cells (blue), infected cells (darker blue), cells
undergoing apoptosis (black), cells undergoing viral load induced pyroptosis
(orange) and cells undergoing bystander induced pyroptosis (red) are shown in
the cell maps, and the corresponding cell counts at each time-point are
displayed in the time-course plots. Note, here $\Pi=1$.  Figure 9: Cell map
progressions and averaged time-course data for three situations, investigating
the probability, $\Pi$, that a cell undergoes pyroptosis, instead of
apoptosis, in response to viral load. We test three values: (A) $\Pi=0$, (B)
$\Pi=0.5$ and (C) $\Pi=1$. Viable cells (blue), infected cells (darker blue),
cells undergoing apoptosis (black), cells undergoing viral induced pyroptosis
(orange) and cells undergoing bystander induced pyroptosis (red) are shown in
the cell maps, and the corresponding cell counts at each time-point are
displayed in the time-course plots. Note, here
$R_{\text{IL-1}\beta}=R_{\text{est}}$, the baseline value.

##  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 [26,
24, 11, 22, 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 trial111Chinese Clinical Trial Registry, registration
number: ChiCDrug2000030002. [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.

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 I 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 [21, 13] and cancer
[26, 23, 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/.

## References

  * [1] P. Zhou, X. Yang, X. Wang, B. Hu, L. Zhang, W. Zhang, H. Si, Y. Zhu, B. Li, C. Huang, et al,  A pneumonia outbreak associated with a new coronavirus of probable bat origin,  Nature, 579 (2020), 270–273.
  * [2] J.K.Y. Yap, M. Moriyama, and A. Iwasaki,  Inflammasomes and pyroptosis as therapeutic targets for COVID-19,  J. Immunol., 205 (2020), 307–312.
  * [3] M.Z. Tay, C.M. Poh, L. Rénia, P.A. MacAry, and L.F.P. Ng,  The trinity of COVID-19: Immunity, inflammation and intervention,  Nat. Rev. Immunol., 579 (2020), 363–374.
  * [4] M. Soy, G. Keser, P. Atagündüz, F. Tabak, I. Atagündüz, and S. Kayhan,  Cytokine storm in COVID-19: Pathogenesis and overview of anti-inflammatory agents used in treatment,  Clin. Rheumatol., 39 (2020), 2085–2094.
  * [5] D. Tang, P. Comish, and R. Kang,  The hallmarks of COVID-19 disease,  PLoS Pathog., 16 (2020), e1008536.
  * [6] A. Shah,  Novel coronavirus-induced NLRP3 inflammasome activation: A potential drug target in the treatment of COVID-19,  Front. Immunol., 11 (2020), 1021.
  * [7] M.Z. Ratajczak and M. Kucia,  SARS-CoV-2 infection and overactivation of NLRP3 inflammasome as a trigger of cytokine storm and risk factor for damage of hematopoietic stem cells,  Leukemia, 4 (2020), 1726-–1729.
  * [8] Y. Fu, Y. Cheng, and Y. Wu,  Understanding SARS-CoV-2-mediated inflammatory responses: From mechanisms to potential therapeutic tools,  Virol. Sin., 35 (2020), 266-–271.
  * [9] S. Nagashima, M.C. Mendes, A.P. Camargo, N.H. Borges, T.M. Godoy, A.F.R. Miggiolaro, F. da Silva Dezidério, C. Machado-Souza, and L. De Noronha,  Endothelial dysfunction and thrombosis in patients with COVID-19-brief report,  Arterioscler. Thromb. Vasc. Biol., 40 (2020), 2404–2407.
  * [10] Y. Jamilloux, T. Henry, A. Belot, S. Viel, M. Fauter, T. El Jammal, T. Walzer, B. François, and P. Sève,  Should we stimulate or suppress immune responses in COVID-19: Cytokine and anti-cytokine interventions,  Autoimmun. Rev., 19 (2020), 102567.
  * [11] C.Y. Taabazuing, M.C. Okondo, and D.A. Bachovchin,  Pyroptosis and apoptosis pathways engage in bidirectional crosstalk in monocytes and macrophages,  Cell Chem. Biol., 24 (2017), 507–514.
  * [12] N. Kelley, D. Jeltema, Y. Duan, and Y. He,  The NLRP3 inflammasome: An overview of mechanisms of activation and regulation,  Int. J. Mol. Sci., 20 (2019), 3328.
  * [13] T. Bergsbaken, S.L. Fink, and B.T. Cookson,  Pyroptosis: Host cell death and inflammation,  Nat. Rev. Microbiol., 7 (2009), 99–109.
  * [14] C.A. Dinarello,  Immunological and inflammatory functions of the interleukin-1 family,  Annu. Rev. Immunol., 27 (2009), 519–550.
  * [15] Y. He, H. Hara, and G. Núñez,  Mechanism and regulation of NLRP3 inflammasome activation,  Trends Biochem. Sci., 41 (2016), 1012–1021.
  * [16] Z.B. Zalinger, R. Elliott, and S.R. Weiss,  Role of the inflammasome-related cytokines IL-1 and IL-18 during infection with murine coronavirus,  J. Neurovirol., 23 (2017), 845–854.
  * [17] A. Stutz, D.T. Golenbock, and E. Latz,  Inflammasomes: Too big to miss,  J. Clin. Invest., 119 (2009), 3502–3511.
  * [18] J.J. O’Shea, M. Gadina, R.M. Siegel, and J. Farber,  Cytokines,  in Rheumatology, (2015), 99–112.
  * [19] P. Song, W. Li, J. Xie, Y. Hou, and C. You,  Cytokine storm induced by SARS-CoV-2,  Clin. Chim. Acta, 509 (2020), 280–287.
  * [20] C. Huang, Y. Wang, X. Li, L. Ren, J. Zhao, Y. Hu, L. Zhang, G. Fan, J. Xu, X. Gu, et al,  Clinical features of patients infected with 2019 novel coronavirus in Wuhan, China,  The Lancet, 395 (2020), 497–506.
  * [21] M. Bertinaria, S. Gastaldi, E. Marini, and M. Giorgis,  Development of covalent NLRP3 inflammasome inhibitors: Chemistry and biological activity,  Arch. Biochem. Biophys., 670 (2019), 116–139.
  * [22] K. Tsuchiya,  Inflammasome-associated cell death: Pyroptosis, apoptosis, and physiological implications,  Microbiol. Immunol., 64 (2020), 252–269.
  * [23] S.L. Fink and B.T. Cookson,  Apoptosis, pyroptosis, and necrosis: Mechanistic description of dead and dying eukaryotic cells,  Infect. Immun., 73 (2005), 1907–1916.
  * [24] A.G. Kozloski,  Inflammasome,  Mater. Methods, 10 (2020), 2869.
  * [25] S. Christgen, D.E. Place, and T.D. Kanneganti,  Toward targeting inflammasomes: Insights into their regulation and activation,  Cell Res., 30 (2020), 315–327.
  * [26] Z. Zheng and G. Li,  Mechanisms and therapeutic regulation of pyroptosis in inflammatory diseases and cancer,  Int. J. Mol. Sci., 21 (2020), 1456.
  * [27] M.G. Sanders, M.J. Parsons, A.G. Howard, J. Liu, S.R. Fassio, J.A. Martinez, and L. Bouchier-Hayes,  Single-cell imaging of inflammatory caspase dimerization reveals differential recruitment to inflammasomes,  Cell Death Dis, 6 (2015), e1813.
  * [28] C. Semino, S. Carta, M. Gattorno, R. Sitia, and A. Rubartelli,  Progressive waves of IL-1$\beta$ release by primary human monocytes via sequential activation of vesicular and gasdermin D-mediated secretory pathways,  Cell Death Dis., 9 (2018), 1–14.
  * [29] G. Lopez-Castejon and D. Brough,  Understanding the mechanism of IL-1$\beta$ secretion,  Cytokine Growth Factor Rev., 22 (2011), 189–195.
  * [30] P. Broz, P. Pelegrín, and F. Shao,  The gasdermins, a protein family executing cell death and inflammation,  Nat. Rev. Immunol., 20 (2019), 143-–157.
  * [31] L. DiPeso, D.X. Ji, R.E. Vance, and J.V. Price,  Cell death and cell lysis are separable events during pyroptosis,  Cell Death Disc., 3 (2017), 1–10.
  * [32] D. Brough and N.J. Rothwell,  Caspase-1-dependent processing of pro-interleukin-1$\beta$ is cytosolic and precedes cell death,  J. Cell Sci., 120 (2007), 772–781.
  * [33] K. Schleich and I.N. Lavrik,  Mathematical modeling of apoptosis,  Cell Comm. Signal., 11 (2013), 1–7.
  * [34] S.L. Spencer and P.K. Sorger,  Measuring and modeling apoptosis in single cells,  Cell, 144 (2011), 926–939.
  * [35] W. Wang and T. Zhang,  Caspase-1-mediated pyroptosis of the predominance for driving CD4 T cells death: A nonlocal spatial mathematical model,  Bull. Math. Biol., 80 (2018), 540–582.
  * [36] D. Veltman, T. Laeremans, E. Passante, and H. J. Huber,  Signal transduction analysis of the NLRP3-inflammasome pathway after cellular damage and its paracrine regulation,  J. Theor. Biol., 415 (2017), 125–136.
  * [37] Y. Bozkurt, A. Demir, B. Erman, and A. Gül,  Unified modeling of familial mediterranean fever and cryopyrin associated periodic syndromes,  Comp. Math. Meth. Med., 15 (2015), 893507.
  * [38] F. López-Caamal and H.J. Huber,  Stable IL-$1\beta$-activation in an inflammasome signalling model depends on positive and negative-feedbacks and tight regulation of protein production,  IEEE ACM T. Comput. Bi., 16 (2018), 627–637.
  * [39] WHO,  World Health Organisation Statement on the Pandemic.  https://www.who.int/dg/speeches/detail/who-director-general-s-opening-remarks-at-the-media-briefing-on-covid-19---11-march-2020, 2020\.
  * [40] F. Akpinar, B. Inankur, and J. Yin,  Spatial-temporal patterns of viral amplification and interference initiated by a single infected cell,  J. Virol., 90 (2016), 7552–7566.
  * [41] A. Bankhead, E. Mancini, A.C. Sims, R.S. Baric, S. McWeeney, and P.M.A. Sloot,  A simulation framework to investigate in vitro viral infection dynamics,  J. Comput. Sci., 4 (2013), 127–134.
  * [42] A.L. Bauer, C. Beauchemin, and A.S. Perelson,  Agent-based modeling of host–pathogen systems: The successes and challenges,  Inf. Sci., 179 (2009), 1379–1389.
  * [43] G. Bocharov, A. Meyerhans, N. Bessonov, S. Trofimchuk, and V. Volpert,  Spatiotemporal dynamics of virus infection spreading in tissues,  PLoS One, 11 (2016), e0168576.
  * [44] C. Beauchemin, S. Forrest, and F.T. Koster,  Modeling influenza viral dynamics in tissue,  in International Conference on Artificial Immune Systems (Springer), (2006), 23–36.
  * [45] C. Beauchemin,  Probing the effects of the well-mixed assumption on viral infection dynamics,  J. Theor. Biol., 242 (2006), 464–477.
  * [46] C. Beauchemin, J. Samuel, and J. Tuszynski,  A simple cellular automaton model for influenza A viral infections,  J. Theor. Biol., 232 (2005), 223–234.
  * [47] D. Levin, S. Forrest, S. Banerjee, C. Clay, J. Cannon, M. Moses, and F. Koster,  A spatial model of the efficiency of T cell search in the influenza-infected lung,  J. Theor. Biol., 398 (2016), 52–63.
  * [48] N. Fachada, V.V. Lopes, and A. Rosa,  Simulating antigenic drift and shift in influenza A,  in Proceedings of the 2009 ACM symposium on Applied Computing, (2009), 2093–2100.
  * [49] A.L. Jenner, F. Frascoli, A.C.F. Coster, and P.S. Kim,  Enhancing oncolytic virotherapy: Observations from a Voronoi cell-based model,  J. Theor. Biol., 485 (2020), 110052.
  * [50] D. Wodarz, A. Hofacre, J.W. Lau, Z. Sun, H. Fan, and N.L. Komarova,  Complex spatial dynamics of oncolytic viruses in vitro: Mathematical and experimental approaches,  PLoS Comput. Biol., 8 (2012), e1002547.
  * [51] G. An,  Introduction of an agent-based multi-scale modular architecture for dynamic knowledge representation of acute inflammation,  Theor. Biol. Med., 5 (2008), 11.
  * [52] R.C. Cockrell and G. An,  Examining the controllability of sepsis using genetic algorithms on an agent-based model of systemic inflammation,  PLoS Comput. Biol., 14 (2018), e1005876.
  * [53] F. Fatehi, R.J. Bingham, E.C. Dykeman, P.G. Stockley, and R. Twarock,  Comparing antiviral strategies against COVID-19 via multi-scale within host modelling,  preprint, arXiv: 2010.08957.
  * [54] T.J. Sego, J.O. Aponte-Serrano, J.F. Gianlupi, S.R. Heaps, K. Breithaupt, L. Brusch, J. Crawshaw, J.M. Osborne, E.M. Quardokus, R.K. Plemper, et al,  A modular framework for multiscale, multicellular, spatiotemporal modeling of acute primary viral infection and immune response in epithelial tissues and its application to drug therapy timing and effectiveness,  PLoS Comput. Biol., 16 (2020), e1008451.
  * [55] Y. Wang, G. An, A. Becker, C. Cockrell, N. Collier, M. Craig, C.L. Davis, J. Faeder, A.N.F. Versypt, J.F. Gianlupi, et al,  Rapid community-driven development of a SARS-CoV-2 tissue simulator,  preprint, bioR$\chi$iv:2020.04.02.019075.
  * [56] A.V. Bagaev, A.Y. Garaeva, E.S. Lebedeva, A.V. Pichugin, R.I. Ataullakhanov, and F.I. Ataullakhanov,  Elevated pre-activation basal level of nuclear NF-$\kappa$B in native macrophages accelerates LPS-induced translocation of cytosolic NF-$\kappa$B into the cell nucleus,  Sci. Reports, 9 (2019), 1–16.
  * [57] C. Zhang, C. Zhao, X. Chen, R. Tao, S. Wang, G. Meng, X. Liu, C. Shao, and X. Su,  Induction of ASC pyroptosis requires gasdermin D or caspase-1/11-dependent mediators and IFN$\beta$ from pyroptotic macrophages,  Cell Death Dis., 11 (2020), 470.
  * [58] Y. Huang, H. Jiang, and Y. Chen, et al,  Tranilast directly targets NLRP3 to treat inflammasome-driven diseases,  EMBO Mol. Med., 10 (2018), e8689.
  * [59] M.S. Salahudeen and P.S. Nishtala,  An overview of pharmacodynamic modelling, ligand-binding approach and its application in clinical practice,  Saudi Pharm. J., 25 (2017), 165–175.
  * [60] N.M. de Vasconcelos, N. Van Opdenbosch, H. Van Gorp, E. Parthoens, and M. Lamkanfi,  Single-cell analysis of pyroptosis dynamics reveals conserved GSDMD-mediated subcellular events that precede plasma membrane rupture,  Cell Death Diff., 26 (2019), 146–161.
  * [61] S. Han, T.B. Lear, J.A. Jerome, S. Rajbhandari, C.A. Snavely, D.L. Gulick, K.F. Gibson, C. Zou, B.B. Chen, and R.K. Mallampalli,  Lipopolysaccharide primes the NALP3 inflammasome by inhibiting its ubiquitination and degradation mediated by the SCFFBXL2 E3 ligase,  J. Bio. Chem., 290 (2015), 18124–18133.
  * [62] J. Cheng, A.L. Waite, E.R., Tkaczyk, K. Ke, N. Richards, A.J. Hunt and D.L. Gumucio,  Kinetic properties of ASC protein aggregation in epithelial cells,  J. Cell. Physio., 222 (2010), 738–747.
  * [63] J. Ruland,  Inflammasome: Putting the pieces together,  Cell, 156 (2014), 1127–1129.
  * [64] J. Chai and Y. Shi,  Apoptosome and inflammasome: Conserved machineries for caspase activation,  Nat. Sci. Rev., 1 (2014), 101–118.
  * [65] A. Iliev, N. Kyurkchiev, and S. Markov,  On the approximation of the step function by some sigmoid functions,  Math. Comput. Simul., 133 (2017), 223–234.
  * [66] MATLAB,  version 1.8.0_202 (R2019n).  The MathWorks Inc., Natick, Massachusetts, 2019.
  * [67] M.A. Moors and S.B. Mizel,  Proteasome-mediated regulation of interleukin-1$\beta$ turnover and export in human monocytes,  J. Leukocyte Biol., 68 (2000), 131–136.
  * [68] F. Martín-Sánchez, C. Diamond, M. Zeitler, A.I. Gomez, A. Baroja-Mazo, J. Bagnall, D. Spiller, M. White, M.J.D. Daniels, A. Mortellaro, et al,  Inflammasome-dependent IL-1$\beta$ release depends upon membrane permeabilisation,  Cell Death Diff., 23 (2016), 1219–1231.
  * [69] G. Qian and A. Mahdi,  Sensitivity analysis methods in the biomedical sciences,  Math Biosci., 323 (2020), 108306.
  * [70] S. Hamis, S. Stratiev and G.G. Powathil,  Uncertainty and sensitivity analyses methods for agent-based mathematical models: An introductory review,  in _The Physics of Cancer: Research Advances_ (ed. Bernard Gerstman), Singapore: World Scientific Publishing, (2021).
  * [71] R.P. Dickinson and R.J. Gelinas,  Sensitivity analysis of ordinary differential equation systems - A direct method,  J. Comput. Phys., 21 (1976), 123–143.
  * [72] A. Ghaffarizadeh, R. Heiland, S.H. Friedman, S.M. Mumenthaler and P. Macklin,  PhysiCell: An open source physics-based cell simulator for 3-D multicellular systems,  PLoS Comput. Biol., 14 (2018), e1005991.
  * [73] K. Lipinska, K.E. Malone, M. Moerland, and C. Kluft,  Applying caspase-1 inhibitors for inflammasome assays in human whole blood,  J. Immunol. Meth., 411 (2014), 66–69.
  * [74] K. Schroder and J. Tschopp,  The inflammasomes,  Cell, 140 (2010), 821–832.
  * [75] W. He, H. Wan, L. Hu, P. Chen, X. Wang, Z. Huang, Z. Yang, C. Zhong, and J. Han,  Gasdermin D is an executor of pyroptosis and required for interleukin-1$\beta$ secretion,  Cell Res., 25 (2015), 1285–1298.
  * [76] S. Lee, M. Hirohama, M. Noguchi, K. Nagata and A. Kawaguchi,  Influenza A virus infection triggers pyroptosis and apoptosis of respiratory epithelial cells through the type I interferon signaling pathway in a mutually exclusive manner,  J. Virol., 92 (2018), e00396-18.
  * [77] M.P. Lythgoe and P. Middleton,  Ongoing clinical trials for the management of the COVID-19 pandemic,  Trends Pharmacol. Sci., 41 (2020), 363–382.
  * [78] K. Sharun, R. Tiwari, J. Dhama, and K. Dhama,  Dexamethasone to combat cytokine storm in COVID-19: Clinical trials and preliminary evidence,  Int. J. Surg., 82 (2020), 179–181.
  * [79] RECOVERY Collaborative Group,  Dexamethasone in hospitalized patients with covid-19: Preliminary report,  N. Engl. J. Med., 384 (2020), 693–704.
  * [80] T. Rhen and J.A. Cidlowski,  Antiinflammatory action of glucocorticoids: New mechanisms for old drugs,  N. Engl. J. Med., 353 (2005), 1711–1723.
  * [81] M. Cornut, E. Bourdonnay and H. Thomas,  Transcriptional Regulation of Inflammasomes,  Int. J. Mol. Sci., 21 (2020), 8087.
  * [82] S. Christgen, D.E. Place and T..D. Kanneganti,  Toward targeting inflammasomes: Insights into their regulation and activation,  Cell Res., 30 (2020), 315–327.
  * [83] F. Chen, G. Jiang, H. Liu, Z. Li, Y. Pei, H. Wang, H. Pan, H. Cui, J. Long, J. Wang, et al,  Melatonin alleviates intervertebral disc degeneration by disrupting the IL-1$\beta$/NF-$\kappa$B-NLRP3 inflammasome positive feedback loop,  Bone Res., 8 (2020), 1–13.
  * [84] H. Wu, C. Zhao, Q. Xie, J. Xu and G. Fei,  TLR2-melatonin feedback loop regulates the activation of NLRP3 inflammasome in murine allergic airway inflammation,  Front. Immunol., 11 (2020), 172.
  * [85] B.B. Mishra, V.A.K. Rathinam, G.W. Martens, A.J. Martinot, H. Kornfeld, K.A. Fitzgerald, and C.M. Sassetti,  Nitric oxide controls the immunopathology of tuberculosis by inhibiting NLRP3 inflammasome–dependent processing of IL-1$\beta$,  Nat. Immunol., 14 (2013), 52–60.
  * [86] G. Guarda, M. Braun, F. Staehli, A. Tardivel, C. Mattmann, I. Förster, M. Farlik, T. Decker, R.A. Du Pasquier, P. Romero, et al,  Type I interferon inhibits interleukin-1 production and inflammasome activation,  Immunity, 34 (2011), 213–223.
  * [87] K.G. Lokugamage, A. Hage, M. de Vries, A.M. Valero-JImenez, C. Schindewolf, M. Dittmann, R. Rajsbaum, and V.D. Menachery,  Type I interferon susceptibility distinguishes SARS-CoV-2 from SARS-CoV,  J. Virol., 94 (2020), e01410-20.
  * [88] C. Juliana, T. Fernandes-Alnemri, S. Kang, A. Farias, F. Qin, and E.S. Alnemri,  Non-transcriptional priming and deubiquitination regulate NLRP3 inflammasome activation,  J. Biol. Chem., 287 (2012), 36617–36622.
  * [89] J.S. Lolkema and D. Slotboom,  The Hill analysis and co-ion–driven transporter kinetics,  J. Gen. Physio., 145 (2015), 565–574.

##  Appendix 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

##  Appendix 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.

Label | Symbol | Component Description
---|---|---
$y_{1}$ | $\text{[NF-$\kappa$B\textsubscript{n}]}(t)$ | concentration of nuclear NF-$\kappa$B
$y_{2}$ | $\text{[NLRP3\textsubscript{i}]}(t)$ | concentration of inactive NLRP3 protein
$y_{3}$ | $\text{[NLRP3\textsubscript{a}]}(t)$ | concentration of active NLRP3 protein
$y_{4}$ | $\text{[NLRP3\textsubscript{o}]}(t)$ | concentration of oligomerised NLRP3 protein
$y_{5}$ | $[\text{ASC\textsubscript{b}}](t)$ | concentration of bound ASC protein
$y_{6}$ | $\text{[C1]}(t)$ | concentration of activated caspase-1
$y_{7}$ | $\text{[GSDMD-N]}(t)$ | concentration of activated GSDMD-N
$y_{8}$ | $[\text{pro-IL-1$\beta$}](t)$ | concentration of pro-IL-1$\beta$
$y_{9}$ | $\text{[IL-1$\beta$\textsubscript{c}]}(t)$ | concentration of active IL-1$\beta$ in the cytoplasm
$y_{10}$ | $\text{[IL-1$\beta$\textsubscript{e}]}(t)$ | concentration of active IL-1$\beta$ external
$y_{11}$ | $[\text{IL-18}_{\text{c}}](t)$ | concentration of active IL-18 in the cytoplasm
$y_{12}$ | $[\text{IL-18}_{\text{e}}](t)$ | concentration of active IL-18 external
$y_{13}$ | $[\text{Drug}](t)$ | concentration of free anti-inflammatory drug
$y_{14}$ | $[\text{Drug}\cdot\text{NLRP3}_{\text{a}}](t)$ | concentration of anti-inflammatory drug - NLRP3${}_{\text{a}}$ complex
$y_{15}$ | $V(t)$ | cell volume
Table B.1: The symbols and descriptions of the components considered in the
model. The bracket notation [ ] denotes compound concentration. Implicit time-
dependence is denoted by $(t)$.

###  B.1 Signal 1.

The initial signal, $S_{1}$, as defined in Equation (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,

$\alpha\frac{[A]^{\gamma}}{A_{50}^{\gamma}+[A]^{\gamma}}.$ (33)

Here, A50 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,

$A\xrightarrow{k}B,$ (34)

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,

$A\xrightarrow{\delta}\emptyset,$ (35)

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,

$A+B\xrightarrow{k}C+B.$ (36)

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,

$\text{[NLRP3\textsubscript{a}]}+\text{[NLRP3\textsubscript{m-1}]}\xrightarrow{k_{2}}\text{[NLRP3\textsubscript{m}]}\quad\mbox{for}\
\text{m}=2,\ldots,\text{p},$

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

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

$\displaystyle\frac{d\text{[NLRP3\textsubscript{m}]}}{dt}$ $\displaystyle=$
$\displaystyle
k_{2}\text{[NLRP3\textsubscript{a}]}\text{[NLRP3\textsubscript{m-1}]}-k_{2}\text{[NLRP3\textsubscript{a}]}\text{[NLRP3\textsubscript{m}]}$
$\displaystyle=$ $\displaystyle
k_{2}\text{[NLRP3\textsubscript{a}]}\left(\text{[NLRP3\textsubscript{m-1}]}-\text{[NLRP3\textsubscript{m}]}\right),$

and

$\frac{d\text{[NLRP3\textsubscript{p}]}}{dt}=k_{2}\text{[NLRP3\textsubscript{a}]}\text{[NLRP3\textsubscript{p-1}]}.$

The total level of bound or oligomerised protein is,

$\text{[NLRP3\textsubscript{o}]}=\sum_{\text{m}=2}^{\text{p}}\text{[NLRP3\textsubscript{m}]}.$

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

$\displaystyle\frac{d\text{[NLRP3\textsubscript{o}]}}{dt}=\sum_{\text{m}=2}^{\text{p}}\frac{d\text{[NLRP3\textsubscript{m}]}}{dt}$
$\displaystyle=k_{2}\text{[NLRP3\textsubscript{a}]}\left(\text{[NLRP3\textsubscript{a}]}-\text{[NLRP3\textsubscript{2}]}+\sum_{\text{m}=3}^{\text{p}-1}\left(\text{[NLRP3\textsubscript{m-1}]}-\text{[NLRP3\textsubscript{m}]}\right)+\text{[NLRP3\textsubscript{p-1}]}\right).$

This simplifies to,

$\frac{d\text{[NLRP3\textsubscript{o}]}}{dt}=k_{2}\text{[NLRP3\textsubscript{a}]}^{2},$

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
NLRP3a, as given in Equation (27). Also included in this modified model are
the dynamics of the anti-inflammatory drug, Drug, and the complex Drug $\cdot$
NLRP3a, as listed in Equations (28) and (29). Here, we consider a binding
reaction of the form,

$\text{Drug +
NLRP3\textsubscript{a}}\underset{k_{\text{-D}}}{\stackrel{{\scriptstyle
k_{\text{+D}}}}{{\rightleftharpoons}}}\text{Drug}\cdot\text{NLRP3\textsubscript{a}}.$
(37)

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
NLRP3o, in an effort to minimise the model.

##  Appendix C Model parameter values

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

Label | Symbol | Description | Value
---|---|---|---
| $a$ | sigmoid constant for step-function approximation | 1 (a.u.)
$\theta_{1}$ | $\alpha_{1}$ | transcription rate of NLRP3i | 0.07 (a.u.) min-1
$\theta_{2}$ | $\alpha_{2}$ | cleavage rate of GDSMD | 0.1 min-1
$\theta_{3}$ | $\alpha_{3}$ | transcription rate of pro-IL-1$\beta$ | 0.06 (a.u.) min-1
$\theta_{4}$ | $\alpha_{4}$ | cleavage rate of pro-IL-1$\beta$ | 1 min-1
$\theta_{5}$ | $\alpha_{5}$ | cleavage rate of pro-IL-18 | 1 min-1
| $b$ | sigmoid constant for step-function approximation | 2 (a.u.)
| $c$ | sigmoid constant for step-function approximation | 1000
$\theta_{6}$ | C150 | half-max value of caspase-1 required for cleavage | 0.3 a.u.
$\theta_{7}$ | $\delta_{1}$ | decay rate of NLRP3i and NLRP3a | 0.002 min-1
$\theta_{8}$ | $\delta_{2}$ | decay rate of pro-IL-1$\beta$ and IL-1$\beta\textsubscript{c}$ | 0.004 min-1
$\theta_{9}$ | $h$ | maximum height elevation of the nuclear NF-$\kappa$B peak | 0.55 a.u.
$\theta_{10}$ | $\gamma_{\text{C1}}$ | Hill function coefficient of caspase-1 mediated cleavage | 2
$\theta_{11}$ | $\gamma_{\text{NF}}$ | Hill function coefficient of nuclear NF-$\kappa$B mediated transcription | 2
$\theta_{12}$ | $k_{1}$ | activation rate of NLRP3i | 0.7 min-1
$\theta_{13}$ | $k_{2}$ | oligomerisation rate of NLRP3a | 1 (a.u.)-1min-1
$\theta_{14}$ | $k_{3}$ | rate for NLRP3a \- ASCf interaction | 0.04 (a.u.)-1 min-1
$\theta_{15}$ | $k_{4}$ | rate for ASCb \- (pro-caspase-1) interactions | 0.03 (a.u.)-1 min-1
$\theta_{16}$ | $k_{5}$ | transport rate of IL-1$\beta$c out of cell | 1 min-1
$\theta_{17}$ | $k_{6}$ | transport rate of IL-18c out of cell | 1 min-1
$\theta_{18}$ | $k_{7}$ | rate at which cell volume increases | 0.2 min-1
| $k_{+\text{D}}$ | forward rate constant for the drug-NLRP3a reaction | 0.005 (a.u.)-1 min-1
| $k_{-\text{D}}$ | reverse rate constant for the drug-NLRP3a reaction | 0.00005 min-1
| $n$ | level of NLRP3${}_{\text{o}}$ required to form the inflammasome base | 1 a.u.
$\theta_{19}$ | NF50 | half-max value of nuclear NF-$\kappa$B required for transcription | 0.3 a.u.
$\theta_{20}$ | $s$ | skewness of the nuclear NF-$\kappa$B peak | 0.8
$\theta_{21}$ | $\tau$ | time when the nuclear NF-$\kappa$B peak occurs | 10 min
| $V_{c}$ | critical volume at which cell ruptures | 1.5 a.u.

Table C.1: Parameters values used within numerical simulations of the model.
The justification for these values is given in Appendix C. The parameters
investigated through sensitivity analysis are assigned a label, $\theta_{n}$.
Here ‘a.u.’ refers to arbitrary units of concentration or volume.

###  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\mbox{m}\mbox{ min}^{-1}$. Therefore we
set,

$\delta_{1}=0.002\ \text{min}^{-1},$

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\mbox{m}\ \mbox{min}^{-1}$. Therefore we set,

$\delta_{2}=0.004\ \text{min}^{-1},$

as we assume measure concentrations as arbitrary units.

Figure C.1: We construct a timeline for the pyroptosis pathway using data from
Bagaev et al. [56], de Vasconcelos et al. [60] and Martín-Sánchez et al. [68].
At 0 minutes (start time), Signal 1 ($S_{1}$) turns on. Thereafter, at 10
minutes, nuclear NF-$\kappa$B levels peak, after which the cytoplasmic-to-
nuclear translocation of NF-$\kappa$B ceases in the mathematical model. 77
minutes into the timeline, enough NLRP3 has been transcribed, synthesised and
activated to form the inflammasome base. Then, at 107 minutes, the GSDMD-N
induced pores start forming, allowing the outflux of IL-1$\beta$ and IL-18, as
well as the influx of extracellular water which causes the cell volume to
start increasing. Finally, at 120 minutes (the end time) the cell membrane
completely ruptures. At this point, around 90% of the pro-IL-1$\beta$ and pro-
IL-18 will have been cleaved and subsequently released to the extracellular
environment.

###  C.2 The timeline of pyroptosis.

We use data from three experimental studies [60, 56, 68] to estimate a
timeline for events regarded in our mathematical model. The values of the
model parameters $\alpha_{1}-\alpha_{5},\ \text{C1}_{50},\
\gamma_{\text{C1}},\ \gamma_{\text{NF}},\ k_{1}-k_{7},\ \text{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 Ca2+ 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 Ca2+
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 NLRP3o 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_{-\text{D}}/k_{+\text{D}}=0.01\text{ a.u.}$,
so that pyroptosis is inhibited for drug concentrations over 0.75 a.u..

##  Appendix 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.

##  Appendix 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 (31) with functions (32)
using the Backward Euler method to implement it in to the PhysiCell framework.
Note, that we do not consider the drug in the PhysiCell setting, that is we
let $[\text{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_{\text{est}}$ and the threshold
level of viral RNA required to initiate viral death is denoted by
$R_{\text{RNA}}$. For the results shown in Figures 8 and 9, we set,

$R_{\text{est}}=200\text{ a.u.}\quad\text{and}\quad R_{\text{RNA}}=2\text{
RNA}.$ (38)

The values of the parameters $\Pi$ and $R_{\text{IL-1}\beta}$ are varies, 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
is 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,

$\displaystyle\frac{dV_{F}}{dt}$ $\displaystyle=$ $\displaystyle-r_{F}V_{F},$
(39) $\displaystyle\frac{dV_{CS}}{dt}$ $\displaystyle=$ $\displaystyle-
r_{CS}V_{CS},$ (40) $\displaystyle\frac{dV_{NS}}{dt}$ $\displaystyle=$
$\displaystyle-r_{NS}V_{NS},$ (41)

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,

$r_{F}=0.05\ \text{min}^{-1},\quad r_{CS}=\frac{1}{60}\ \text{min}^{-1},\quad
r_{NS}=\frac{0.35}{60}\ \text{min}^{-1},$ $T_{A}=8.6\ \text{hrs}=516\
\text{mins},\quad V_{m}=2\ \mu\text{m}^{3}.$

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_{\text{RNA}}$.

Figure E.1:  Numerical simulations results of the Backward Euler
discretisation of the ODE model described by system (31) with functions (32).
The parameters are displayed in Table C.1 and here we use a time-step of
$\Delta_{t}=1$ min. We display the concentration of each model component in
arbitrary units (a.u.) over time, as well as cell volume dynamics. Note that
the time progression stops once the critical volume $V_{c}$ is reached and the
ultimate cell membrane rupture occurs. Figure E.2:  The averaged time-course
data for all cases from Figure 8 and Figure 9 with viable cells omitted. Cells
undergoing apoptosis (black), cells undergoing viral induced pyroptosis
(orange) and cells undergoing bystander induced pyroptosis (red) are shown in
the plots. The average over 10 runs is plotted as an opaque line, where the
transparent colour represents the standard deviation.