Innate Immunity Plays a Key Role in Controlling Viral Load in COVID-19: Mechanistic Insights from a Whole-Body Infection Dynamics Model.

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a pathogen of immense public health concern. Efforts to control the disease have only proven mildly successful, and the disease will likely continue to cause excessive fatalities until effective preventative measures (such as a vaccine) are developed. To develop disease management strategies, a better understanding of SARS-CoV-2 pathogenesis and population susceptibility to infection are needed. To this end, mathematical modeling can provide a robust in silico tool to understand COVID-19 pathophysiology and the in vivo dynamics of SARS-CoV-2. Guided by ACE2-tropism (ACE2 receptor dependency for infection) of the virus and by incorporating cellular-scale viral dynamics and innate and adaptive immune responses, we have developed a multiscale mechanistic model for simulating the time-dependent evolution of viral load distribution in susceptible organs of the body (respiratory tract, gut, liver, spleen, heart, kidneys, and brain). Following parameter quantification with in vivo and clinical data, we used the model to simulate viral load progression in a virtual patient with varying degrees of compromised immune status. Further, we ranked model parameters through sensitivity analysis for their significance in governing clearance of viral load to understand the effects of physiological factors and underlying conditions on viral load dynamics. Antiviral drug therapy, interferon therapy, and their combination were simulated to study the effects on viral load kinetics of SARS-CoV-2. The model revealed the dominant role of innate immunity (specifically interferons and resident macrophages) in controlling viral load, and the importance of timing when initiating therapy after infection.

In January 2020, severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) was identified as the infectious agent causing an outbreak of viral pneumonia in Wuhan, China. It was soon established that droplet-based human to human transmission was occurring, and on March 11, 2020, the World Health Organization characterized coronavirus disease 2019 (COVID-19) as a pandemic. As of December 25, 2020, COVID-19 has infected more than 79 million people and caused more than 1.74 million deaths worldwide.[1] A pandemic-scale outbreak creates tremendous socioeconomic burden due to thwarted productivity, a spike in healthcare expenses, and irreparable loss of human lives.[2,3] Furthermore, implementation of social and physical isolation measures has caused many countries to declare states of emergency and lockdowns with border closures.

SARS-CoV-2 is the seventh identified human coronavirus and the third novel one to emerge in the last 20 years. It has a single-stranded positive sense RNA genome of about 30,000 nucleotides that encodes ∼27 proteins including 4 structural proteins. A surface-expressed spike protein mediates receptor binding and membrane fusion with host cells, and the virus interacts with the angiotensin converting enzyme 2 (ACE2) receptor to gain entry into cells.[4] ACE2 mRNA is present in almost all human organs, but the receptor is particularly highly expressed on the surface of lung alveolar epithelial cells and enterocytes of the small intestine, thereby allowing a preferential accumulation of the virus in these organs.[5] The incubation period of SARS-CoV-2 ranges from about 3–17 days, and COVID-19 diagnosis cannot be made based on symptoms alone, as most are nonspecific and may be confused for more common ailments. The more serious sequelae of infection includes acute respiratory distress syndrome (ARDS) and sepsis caused by the cytokine storm from the immune response to infection, which is believed to be the leading cause of mortality in COVID-19 patients.[6] Screening for COVID-19 is done via nucleic acid testing by RT-PCR (specimens from both upper and lower respiratory tracts) and pulmonary CT scans. The viral load in naso- or oro-pharyngeal swabs is the key clinical biomarker of COVID-19 and also the key clinical end point of pharmacological intervention.

Although antiviral, antibody, and immunomodulatory drugs are being used for treatment of various aspects of the infection, there are still few effective therapeutics for COVID-19 to date. To explore novel and effective therapeutic targets, we require a better understanding of the pathogenesis of COVID-19, particularly of virus–host interactions.[7] This will also enable more efficient disease management strategies, such as deriving prognostic information from viral load kinetics, and quantification of the effects of the immune system in controlling the disease. With limited studies on the in vivo dynamics of SARS-CoV-2, a mathematical modeling approach can be an excellent, complementary tool for investigating viral–host interactions and host immune response in order to better understand COVID-19 progression and evaluate treatment strategies. Indeed, the application of mathematical modeling and quantitative methods has been instrumental in our understanding of viral–host interactions of various viruses, including influenza, HIV, HBV, and HCV.[8] The kinetic models for the above viruses have been developed for various spatial scales, including molecular, cellular, multicellular, organ, and organism. By analyzing viral load kinetics, these models have deepened our understanding of the fundamentals of virus–host interaction dynamics, innate and acquired immune response, drug mechanisms of action, and drug resistance.[9−13]

While the fundamental principles governing different viral infections are similar among most viral species, the kinetics of the underlying mechanisms may vary based on the virus type. Researchers are already using mathematical models to understand the outbreak of COVID-19 in order to guide the efforts of governments worldwide in containing the spread of infection. While most of the models developed to date have focused on the epidemiological aspects of COVID-19 to understand the interhuman transmission dynamics of SARS-CoV-2,[14−18] there are a few studies that have investigated its virus–host interactions and pathogenesis. For example, Goyal et al. developed a mathematical model to predict the therapeutic outcomes of various COVID-19 treatment strategies.[19] Their model is based on target cell-limited viral dynamics[20] and incorporates the immune response to infection in order to predict viral load dynamics in patients pre- and post-treatment with various antiviral drugs. This model was used to project viral dynamics under hypothetical clinical scenarios involving drugs with varying potencies, different treatment timings postinfection, and levels of drug resistance, and the results of this study suggest the application of potent antiviral drugs prior to the peak viral load stage, i.e., in the presymptomatic stage, as an effective means of controlling infection in the body. Further, Wang et al. developed a prototype multiscale model to simulate SARS-CoV-2 dynamics at the tissue scale,[7] wherein an agent-based modeling approach was used to simulate intracellular viral replication and spread of infection to neighboring cells. To unravel the mechanistic underpinnings of clinical phenotypes of COVID-19, Sahoo et al. developed a mechanistic model that studies the intercellular interactions between infected cells and immune cells.[21] Also, Ke et al. developed a model to quantify the early dynamics of SARS-CoV-2 infection in the upper and lower respiratory tracts and used the model to predict infectiousness and disease severity based on viral load dynamics and immune response to infection.[22] Although also a target cell-limited model, by only including upper and lower respiratory tract compartments, this model omits key biological mechanisms involved in the complete immune response and is thus unable to provide deeper insights into the system-wide dynamics and interplay of disease response.

In order to improve upon the existing models, we have developed a multiscale semimechanistic model of viral dynamics, which, in addition to capturing virus-host interactions locally, is also capable of simulating the whole-body dynamics of SARS-CoV-2 infection and is thereby capable of providing insights into disease pathophysiology and the typical and atypical presentations of COVID-19. Importantly, using our modeling platform, we can identify treatment strategies for effective viral load suppression under various clinically relevant scenarios. We note that while the modeling platform is developed for SARS-CoV-2, we also expect it to be applicable to other viruses that have shared similarities in mechanisms of infection and physical dimensions.

Results and Discussion

Model Development, Calibration, and Baseline Solution

We have developed a semimechanistic mathematical model to simulate the whole-body biodistribution kinetics of SARS-CoV-2 following infection through the nasal route (Figure ; Methods: Model development). The model was formulated as a system of ordinary differential equations (eqs –40) that describe cellular-scale viral dynamics, whole-body transport and excretion of viruses, and innate and adaptive immune response to predict the viral load kinetics of SARS-CoV-2 in the respiratory tract, plasma, and other organs of the body. SARS-CoV-2 exhibits ACE2-tropism;[23] therefore, the organs included in the model were chosen based on the presence of ACE2 receptor expressing cells.[24−26] Specifically, the key processes described by the model include infection of ACE2-expressing susceptible cells by SARS-CoV-2 (also referred to as target cells), production of new virions by infected cells, death of infected cells due to cytopathic effects, transport of virions from the site of infection to other organs of the body, hepatobiliary excretion of the virions, and key processes in the innate and adaptive immune response against the virus and infected cells to clear the infection. Note that in the absence of a thorough understanding of the mechanistic underpinnings of viral shedding in the feces,[27] and a growing evidence of liver damage in COVID-19 patients,[28,29] we used bile production as the rate-limiting step for fecal excretion of the virus. Due to unavailability of in vivo or clinical data for viral excretion, it was not feasible to accurately estimate the excretion rate parameter; hence, we mechanistically modeled the process of fecal excretion by using the bile production rate as a surrogate for fecal excretion rate.

Figure 1

Model schematic showing system interactions. Connectivity diagram between compartments indicating viral transport mechanisms, cell populations, immune system agents, and their interactions. Estimated characteristic times of the various transport processes are given in parentheses alongside red arrows. Notation: (V) virus, (H) healthy cells, (I) infected cells, (IFN) interferon, (Ab) antibody, (CD8*) effector CD8+ cells, and (APC) antigen presenting cells. Solid red arrows indicate transport of virus; solid green arrows indicate transformation of a cell into another type; solid black arrows indicate production of an agent; purple dashed lines indicate interaction between two agents; and solid dark red arrows with a flat head indicate inhibition. MPS stands for mononuclear phagocytic system and comprises liver and spleen.

While some of the parameters of the model were known a priori (Table ), the remaining parameters were estimated through nonlinear regression using published in vivo(30) and clinical[31] data. Specifically, from published experimental data for hamsters,[30] we first calibrated a reduced version of the model (referred to as Reduced model; eqs –23) that comprises all compartments and interferon (IFN)-mediated innate immunity but lacks adaptive immunity (bottom half of Figure ; also see workflow in Figure ). The parameters of the reduced model characterize cellular-scale viral dynamics, IFN-mediated immunity, intercompartment viral transport, and hepatobiliary excretion of the virus from the mononuclear-phagocytic system (MPS). Note that in this model, the MPS compartment is the liver and spleen combined together. The estimated parameters were then used in the complete version of the model (referred to as Full model; eqs –40), which also includes adaptive immunity, to calibrate the remaining parameters using nonlinear regression with clinical data.[31] The models were solved numerically in MATLAB as an initial value problem, using the built-in stiff ODE solver ode15s.

Table 1

List of Model Parameters Known A Priori

Notation	Definition	Units	Value (Hamster)	Value (Human)	ref
Immunity parameters
DCD4	Death rate of activated CD4 cells	d–1		0.0185	(68,69)
DPS	Death rate of short-lived plasma cells	d–1		0.5	(70−72)
DPL	Death rate of long-lived plasma cells	d–1		0.0083	(71,72)
λCD8	Death rate of activated CD8 cells	d–1		0.0139	(73,69,68)
DAPC	Death rate for activated APCs	d–1		0.2	(74−76)
APC	Generic steady state concentration of APCs	cells·mL–1		105.3	(77)
APCL	LRT steady state concentration of APCs	cells·mL–1		106	(77,78)
APCM	MPS steady state concentration of APCs	cells·mL–1		107.65	(77,78)
CD8	Naïve CD8 steady state concentration	cells·mL–1		105	(79)
BC	Naïve B cells steady state concentration	cells·mL–1		105	(80,79)
CD4	Naïve CD4 steady state concentration	cells·mL–1		105.8	(79)
GAPC	Reequilibration rate of APCs	cell·mL·d		0.1	(74,81)
GCD8	Reequilibration rate of naïve CD8 cells	cell·mL·d		0.0045	(68,69)
GBC	Reequilibration rate of naïve B cells	cell·mL·d		1.75	(82)
GCD4	Reequilibration rate of naïve CD4 cells	cell·mL·d		0.0027	(69,68)
GIFN	Baseline IFN production rate	IFN units·mL–1·d–1	38.8	1.37	(30,31)
DIFN	Degradation rate of IFN	d–1	24	24	(83)
Healthy tissue parameters
yU0	Initial concentration of target cells in URT	cells·mL–1	106.9	106.7	(78,84,85)
yL0	Initial concentration of target cells in LRT	cells·mL–1	106.9	106.7	(84,85)
yG0	Initial concentration of target cells in GI	cells·mL–1	107.4	106.9	(86)
yM0	Initial concentration of target cells in MPS	cells·mL–1	106.8	106.6	(87−89)
yH0	Initial concentration of target cells in heart	cells·mL–1	106.5	106.3	(90,91)
yK0	Initial concentration of target cells in kidneys	cells·mL–1	107.1	106.9	(87,92)
yB0	Initial concentration of target cells in brain	cells·mL–1	105.9	105.7	(93,94)
Viral dynamics parameters
vU0	Initial concentration of viral load in URT	virus·mL–1	105.7	Table 3	(30)
Transport parameters
ClAb	Clearance rate of antibodies	d–1		0.04	(11)
E	Hepatobiliary excretion rate	d–1	2.68	0.15	(95)

Figure 2

Modeling workflow. In the first step, the reduced model uses the hamster data to fit parameters corresponding to viral dynamics, transport coefficients, and innate immunity. Subsequently, the more complex full model uses the previously estimated parameters and clinical data to estimate the remaining parameters pertaining to the adaptive and innate immune responses. Lastly, simulations are generated to make predictions, and statistical analyses are performed to extract useful information.

Calibrating Parameters of the Reduced Model

As shown in Figure , the numerical solution of the reduced model for whole-body viral kinetics, IFN kinetics, and target cell population kinetics in hamsters satisfies the initial conditions and is in good agreement with the available in vivo data[30] for viral and IFN kinetics (Pearson correlation coefficient R between experimental data and model fits is >0.97, p < 0.0001, Figure S1a). The corresponding parameter estimates along with their 90% confidence intervals (CI) are given in Table . The Akaike information criterion (AIC) and Bayesian information criterion (BIC) values for the reduced model are 13.5 and 36.7, respectively. Based on findings in the in vivo study by Chan et al.[30] that the adaptive immune response in test animals was not triggered during the first 7 days postinfection, it is reasonable to use the reduced form of the model to estimate the unknown parameters, rather than using the full model at this point.

Figure 3

Whole-body reduced model calibration with in vivo data. Nonlinear least-squares regression of the whole-body reduced model using published in vivo data for hamsters is shown. The compartments under consideration are (a) upper respiratory tract (URT), (b) lower respiratory tract (LRT), (c) gastrointestinal tract (GI), (d) MPS, (e) kidneys, (f) heart, (g) brain, and (h) plasma. Orange solid and orange dashed lines indicate the population of healthy and infected ACE2-expressing target cells in a given compartment, respectively; the blue solid line represents the viral load in the corresponding compartment. The left y-axis represents viral load, and the right y-axis represents cell populations in (a) through (g). (i) Logarithm base 10 of IFN mRNA copies per γ-actin multiplied by 106 is shown. The triangular markers represent experimental data presented as mean ± SD (n = 5 animals per time point). Pearson correlation coefficient between observed and fitted values is R > 0.97 (Figure S1a).

Table 2

List of Model Parameters Estimated from Fitting the Reduced Model to In Vivo Dataa

Notation	Definition	Units	Value (90% CI)
Viral dynamics parameters
I	Infection rate of target cells	mL·virus–1·d–1	4.05 (0–81.1)
Pv	Virus production rate in infected cells	virus·cell–1·d–1	19.03 (0–696)
D1	Cytopathic death rate of infected cells	d–1	0.15 (0–0.49)
Innate immunity parameters
PIFN	IFN production rate from infected cells	mRNA·mL·γ actin–1 cell–1·d–1	3.67 (0–8.1)
ε	IFN efficacy	mRNA–1·γ actin–	2.59 (0–103)
λLMϕ	Viral elimination rate by macrophages in LRT	d–1	1.31 (0–9.3)
λMMϕ	Viral elimination rate by macrophages in MPS	d–1	4.24 (0–19.6)
Transport parameters
C1	URT to LRT viral transport rate	d–1	0.73 (0.16–1.29)
C2	URT to GI viral transport rate	d–1	0.0016 (0–0.12)
A1	Viral pulmonary absorption rate	d–1	0.33 (0–8.6)
A2	Viral intestinal absorption rate	d–1	1.24 (0.16–2.3)
HA	Hepatic artery-mediated viral transport rate	d–1	19.2 (0–626)
HV	Hepatic vein-mediated viral transport rate	d–1	17.8 (0–715)
RA	Renal artery-mediated viral transport rate	d–1	8.5 (0–133)
RV	Renal vein-mediated viral transport rate	d–1	12.07 (0–165)
COA	Coronary artery-mediated viral transport rate	d–1	7.2 (0–167)
COV	Coronary vein-mediated viral transport rate	d–1	9.3 (0–186)
CAA	Carotid artery-mediated viral transport rate	d–1	0.97 (0–58.5)
JV	Jugular vein-mediated viral transport rate	d–1	2.32 (0–69.9)

90% confidence intervals (CI) of the estimated parameters are shown in parentheses. For practical purposes, the negative value of the lower bounds of CIs has been replaced by zero. Note: Characteristic times corresponding to the transport parameters are shown in Figure .

The model solution (Figure ) shows the kinetics of ACE2-expressing target cells (solid orange lines) and their infected counterparts (dashed orange lines) in every compartment. These infected cells can produce new virions that will in turn infect other healthy target cells. Because we are using a target-cell limited modeling assumption,[10,19] the healthy target cells that become infected by the virus are not replaced by new healthy cells, and as shown in Figure , the target cells were observed to deplete within 48 h postinfection. The viral load kinetics (blue curve) is primarily governed by the interplay of new virion production, distribution of the virions between compartments, viral elimination by alveolar and MPS macrophages, hepatobiliary excretion of viruses from the body, cytopathic death of infected cells, and suppression of viral production due to IFN produced by infected cells,[10] which is shown in Figure i. As the infected cell population tapers, the IFN concentration will also decrease to the preinfection baseline value. In our model, infected respiratory tract cells are the source of IFN following infection, the lack of which has been found to be the underlying cause of life-threatening COVID-19 due to uncontrolled viral replication in the absence of IFN regulation.[32−34]

Of note, the plasma compartment (Figure h) of the model does not contain any target cell population and thus its viral load kinetics is only governed by the influx and outflux of viruses from various compartments. However, in the full model, the neutralization of viruses by antibodies will also be considered in the plasma compartment, as discussed in the next section. Plasma flow is the key mechanism of viral transport and systemic spread of infection in the body,[35] but due to lack of established mechanistic underpinnings of these processes, we instead use phenomenological rate constants to characterize viral transport. Based on the estimated characteristic times (1–24 h) of the vascular transport processes (shown in Figure and presented as rates in Table ), it can be inferred that viral transport is permeability-limited and not perfusion-limited, i.e., capillary permeability and vascular surface area govern the rate of extravasation of virions from blood vessels into tissue interstitium to reach the target cells, and thus viral transport is not exclusively governed by the plasma flow rates into the organs. This is consistent with the in vivo behavior of nanomaterials of comparable size[36−44] and is in contrast to the perfusion rate-limited kinetics of smaller lipophilic molecules. The variability in characteristic times of vascular transport can be explained by differences in the permeability of capillary endothelium due to differences in pore sizes of endothelial fenestrae.[45] For instance, the blood–brain barrier seems to resist transport of virus to the brain, thereby leading to an estimated characteristic influx time of 1 day, which is ∼20 times longer than the estimated characteristic time of influx to the MPS (1.25 h), which contains large sinusoidal pores in its microvasculature. Of note, the nonvascular transport processes have relatively longer characteristic times that can be attributed to resistance to transport offered by mucus or degradation caused by pH, among other factors.

Calibrating Parameters of the Full Model

Once the parameters discussed in the previous section were estimated, they were then used in the full model to calibrate the remaining parameters (see Table ) relevant to the innate and adaptive immune system using published clinical data (n = 4 untreated patients with mild symptoms).[31] Due to the uncertainty associated with the duration between time of infection and onset of symptoms (referred to as incubation period), a shifting parameter τ was included in the calibration routine. Numerically, the time points corresponding to the data were shifted τ units of time. As shown in Figure , the model correctly represents the initial conditions of the variables and predicts an incubation period τ of ∼6 days (indicated by the red arrow in Figure a), which is comparable to published literature.[46] Also, assuming the nasal route as the route of infection, the numerical solution for the upper respiratory tract (URT) at time t = 0 suggests exposure to a viral load of ∼107 copies mL–1. The clinical data shows the viral load kinetics in the upper (Figure a) and lower respiratory tract (Figure b), the IFN kinetics (Figure i), the effector CD8+ (CD8*) and activated CD4+ (CD4*) cell population kinetics (Figure k), and the total neutralizing antibody kinetics (Figure l). As shown in Figure , the full model solution fits the data well (Pearson correlation coefficient R > 0.98, p < 0.0001, Figure S1b) and predicts the kinetics of viral load in the remaining compartments by using the viral dynamics and transport parameters estimated from the in vivo data (through calibration of the reduced model). The AIC and BIC values for the full model are −3.04 and 24.5, respectively. The model predicts that the viral load in extrapulmonary organs and plasma persists for ∼17–20 days post onset of symptoms, consistent with published studies[47] and is thus comparable to the duration of viral detection in URT and lower respiratory tract (LRT), therefore suggesting that nasopharyngeal and oropharyngeal swabs can perhaps safely indicate the systemic level of viral load in the body. However, this modeling observation requires further experimental validation.

Table 3

List of Model Parameters Estimated from Fitting the Full Model to Clinical Dataa

Notation	Definition	Units	Value (90% CI)
Innate immunity parameters
PIFN	IFN production rate from infected cells	pg·cell–1·d–1	5.37 (0.7–10)
ε	IFN effectiveness	mL·pg–1	4.57 (0.5–8.6)
Viral dynamics parameters
vU0	Initial concentration of viral load in the URT	virus·mL–1	106.9 (105.8 – 107.9)
τ	Viral incubation time	d	5.77 (1.9–9.6)
Adaptive immunity parameters
TAPC	Transition rate of naïve APCs into APC*	mL·virus–1·d–1	65.2 (0–1.3e+3)
TT-cells	Transition rate of naïve CD8 into CD8* or CD4 into CD4*	mL·cell–1·d–1	0.0062 (0.005–0.007)
TBC	Transition rate of naïve B cells into B*	mL·cell–1·d–1	0.029 (0–0.15)
TBCS	Transition rate of B* into short-lived plasma cells	mL·cell–1·d–1	564.8 (0–2.1e+4)
TBCL	Transition rate of B* into long-lived plasma cells	mL·cell–1·d–1	0.36 (0–16.5)
GPS	Production rate of antibodies from short-lived plasma cells	antibody·cell–1·d–1	0.59 (0–1.35)
GPL	Production rate of antibodies from long-lived plasma cells	antibody·cell–1·d–1	1.34 (0–90.2)
λAb	Antibody loss rate due to viruses	mL·virus–1·d–1	0.0014 (0–0.004)
DAb	Elimination rate of viruses due to antibodies	mL·antibody–1·d–1	0.0085 (0–0.04)
DCD8	Death rate of infected cells due to CD8*	mL·cell–1·d–1	0.0037 (0–0.01)

90% confidence intervals (CI) of the estimated parameters are shown in parentheses. For practical purposes, the negative value of the lower bounds of CIs has been replaced by zero.

Figure 4

Whole-body full model calibration with clinical data. Numerical results for the full model, corresponding to the second step given in the modeling workflow in Figure are shown. The compartments under consideration are (a) URT, (b) LRT, (c) GI, (d) MPS, (e) heart, (f) kidneys, (g) brain, (h) plasma, and (j, k, l) lymphatic compartment. Orange solid and orange dashed lines indicate the population of healthy and infected ACE2-expressing target cells in a given compartment, respectively, while the blue solid line represents the viral load in the corresponding compartment. The left y-axis represents viral load, and the right y-axis represents cell populations in panels (a) through (g). (i) Concentration kinetics of IFN. (j and k) Kinetics of adaptive immune system cells in the lymphatic compartment. (j) Naïve/immature CD4+, CD8+, and B cells. (k) Cells in their activated/effector state along with antibody producing plasma cells. (l) Concentration kinetics of total antibody (IgG, IgM, IgA). The red arrow in panel (a) indicates time of onset of symptoms. A lower bound at 102 copies mL–1 and 104.6 pg mL–1 is imposed to represent the detectable limit of viral load and antibodies, respectively; the dashed blue line indicates the limit of detection. Once the viral load or antibodies go below the detection limit, a vertical line is used to indicate the time of occurrence of this event. The triangular markers represent clinical data presented as mean ± SD (n = 4 patients). The Pearson correlation coefficient between observed and fitted values is R > 0.98 (see Figure S1b).

The model also shows the kinetics of naïve lymphocytes (Figure j) and antibody producing plasma cells (Figure k) in the lymphatic compartment, which is represented as a common compartment for the entire body. Importantly, in close agreement with published literature,[48,49] the model predicts that the systemic concentration of antibodies persists above the detectable limit for >100 days post onset of symptoms, following which it may no longer be detectable (Figure l). This finding suggests the lack of indefinite antibody protection against reinfection[50] and highlights the need for vaccine boosters to achieve long-lasting immunity.[51] We note that the data used in the above calibrations was obtained under conditions where neither the animals nor the patients were given any pharmacological treatment. Hence, the data are appropriate to calibrate the effects of the immune components and other physiological processes in distributing and eliminating the viral load.

While URT and LRT are the preferred sites to detect the presence of SARS-CoV-2, it is important to note, and as is evident from the model predictions, that the viral load in nonpulmonary organs can attain comparable levels and may thus explain the nonrespiratory symptoms observed in some COVID-19 patients.[29,52−56] Following transport of the virus from the respiratory tract to blood or via the gastrointestinal tract to blood, organs that have a significant population of ACE2-expressing cells may become infected by the virus, leading to the extrapulmonary manifestations of COVID-19 that may include symptoms such as diarrhea and impaired renal-, hepatic-, cardiovascular-, or neurological-functions. This suggested relationship between organ viral load and extrapulmonary symptoms demands further experimental investigation and clinical validation.

Individualized Effects of Immune Components on Viral Kinetics

To investigate the effects of the cellular and humoral arms of innate and adaptive immunity on viral load kinetics in the body, we individually switched off these components and simulated the whole-body viral kinetics up to the time when viral load fell below the detectable limit[19] of 102 copies mL–1. This numerical experiment is meant to mimic the effect of compromised immunity due to an underlying condition in a virtual patient undergoing no antiviral treatment. As seen in the viral kinetics shown in Figure , when one or all of the immune components were shut down, the viral concentration in all the compartments was higher than the baseline (dashed dark red line). Further, it can be inferred that IFN is the primary mechanism of controlling viral load in the URT (Figure a) and GI (Figure c), while macrophages (alveolar macrophages, Kupffer cells, and splenic macrophages) play a predominant role in limiting infection in the LRT, MPS, plasma, and other organs, followed by IFNs (Figures b, 5d–5h). Findings in the literature support the above observations as follows. Lack of IFNs can lead to excessive viral production[34] and cause life-threating COVID-19 in patients deficient in functional IFNs due to, for example, the occurrence of loss-of-function mutations in genes governing IFN-mediated immunity[33] or autoantibodies against IFNs.[32] Further, FABP4+ alveolar macrophages were observed to be largely absent in patients with severe COVID-19 but were the predominant macrophage in patients with mild disease,[57] indicating the major role of FABP4+ alveolar macrophages in controlling infection as also shown previously for patients with chronic obstructive pulmonary disease (COPD).[58] We assume a constant supply of macrophages in the lungs and MPS in our model, but a deleterious effect of infection on these immune cells cannot be ruled out, and further experimental evidence is necessary to model the macrophage population kinetics more appropriately.[59] Further, the model reveals that the effect of adaptive immunity (antibodies and CD8* cells) is not significant in controlling infection, but it does not necessarily rule out the therapeutic potential of exogenously administered antibodies or novel cell-based therapies (e.g., T cell therapy).

Figure 5

Effects of different components of the immune system. Comparison between different scenarios where a particular mechanism of immunity has been turned off is shown. The viral kinetics in a given compartment is presented on the y-axis. The compartments under consideration are (a) URT, (b) LRT, (c) GI, (d) MPS, (e) heart, (f) kidneys, (g) brain, and (h) plasma. Six scenarios are considered. (1) Baseline: All immune mechanisms are functional. (2) IFN(−): the innate immune effects of IFN are deactivated. (3) Macrophage(−): alveolar and MPS macrophages are eliminated. (4) Ab(−): no antibodies are produced. (5) CD8*(−): cytotoxic effects of CD8* are suppressed. (6) No immunity: all immunity mechanisms are absent. The flat line at 1012 copies mL–1 represents the upper bound imposed on viral load to showcase only clinically relevant results. A lower bound at 102 copies mL–1 is imposed to represent the detectable limit of viral load in the body, indicated by the horizontal black line; once the viral load goes below the detection limit, a vertical line is used to indicate the time of occurrence of this event.

All the above scenarios abstractly represent real-world underlying conditions (e.g., cancer, diabetes, autoimmune diseases) in patients that lead to varying degrees of immunosuppression, thus highlighting the importance of an individual’s immune status in regulating SARS-CoV-2 kinetics. In the absence of immune responses, the only plausible mechanism that brings the viral load down is hepatobiliary excretion of the virus through the MPS, but our results show that it takes several weeks before the viral load falls below the detectable value of 102 copies mL–1. Note that upon shutting down IFN- and macrophage-mediated immunity or total immunity, the viral load grew beyond clinically observed values in the literature (∼1012 copies mL–1); therefore, we set an upper bound at 1012 copies mL–1 to keep the results clinically meaningful.

Parametric Analysis

To identify the significance of model parameters in affecting viral load kinetics, we performed global sensitivity analysis using the full model (Figure ). The model outputs used to investigate the influence of perturbations to input model parameters were area under the curves of viral load kinetics in URT, LRT, and plasma from 0 to 30 days (i.e., AUC0–30URT, AUC0–30LRT, and AUC0–30Plasma, respectively) and time to reduce the viral load to <102 copies mL–1 in URT, LRT, and plasma (i.e., tclearURT, tclearLRT, and tclearPlasma, respectively).

Figure 6

Global sensitivity analysis. All parameters are perturbed simultaneously from their reference value. (a−f) The ranking of each parameter is quantified as the sensitivity index (SI, y-axis) measured under six criteria. The first three are area under the curves of viral load kinetics (AUC0–30) for the first 30 days in the (a) URT, (b) LRT, and (c) plasma compartment. The next three are the total time (tclear) required for the viral load to fall below the detectable concentration of 102 copies mL–1 in the (d) URT, (e) LRT, and (f) plasma compartment. Each SI is accompanied by an error bar resulting from the sampling distribution. (g) Heat map indicating the ranking of the parameters based on the magnitude of the sensitivity index.

As shown in Figure a–f, multivariate linear regression analysis (MLRA)-based sensitivity indices provide insight into the relative importance of model parameters in governing viral load kinetics. The ranking obtained from MLRA (Figure g) suggests that the total viral load (AUC0–30) in URT, LRT, and plasma (indicator of systemic behavior) is strongly dependent on viral production rate (P), IFN production rate (PIFN), and cytopathic death rate of infected cells (D1). Cytopathic death of infected cells is a property of the system that may be hard to manipulate pharmacologically; therefore, IFN concentration and viral production rate are conceivably the targets best-suited for pharmacological interventions to constrain the viral load in patients.[60] Treatment over a combination of the two parameters, i.e., IFN and viral production rates, can possibly be more efficient than monotherapy in suppressing viral load, but clinical validation of this hypothesis is necessary. While cytopathic death rate and viral production rate play major roles in governing the clearance time of viruses (tclear), IFN production rate is relatively less significant. Further, transport processes play a significant role in governing both total viral load and time to clear the load. Specifically, the transport of virus from URT to LRT (C1) affects viral load in URT and viral transport from blood to MPS via the hepatic artery (H) and MPS to blood via the hepatic vein (H) govern the clearance time from plasma, which should also affect the viral clearance of other organs connected to plasma. Also, the elimination rates of viruses by alveolar macrophages (λMϕ) and MPS macrophages (λMϕ) are important parameters in governing the kinetics of viral load in LRT and plasma, respectively. Finally, the total viral load in LRT is also affected by pulmonary absorption (A1), i.e., by the rate of transport of viruses from LRT to plasma.

Clinical Application of Model

Up to this point, we calibrated the reduced model with in vivo data and used the estimated parameters to calibrate the remaining parameters of the full model with clinical data for patients with mild symptoms. While the full model was able to produce observations consistent with published literature, the validation of viral load kinetics of the extrapulmonary organs needs corroboration from future clinical studies. However, to demonstrate the possible clinical utility of our model, we next sought to examine whether this tool can make predictions that might allow clinicians to design an effective therapy for patients, helping to optimize and personalize their treatment regimen. As a numerical experiment, we tested three treatment scenarios in controlling infection: a hypothetical antiviral agent, interferon therapy, and antiviral agent–interferon combination therapy. Further, the effects of the timing of therapy initiation (t) were also tested, i.e., starting therapy on the day of onset of symptoms (t = tonset) and starting therapy 5 days post onset of symptoms (t = tonset + 5). To quantify treatment effectiveness, we compared the viral load kinetics in three compartments, namely, URT, LRT, and plasma in simulations spanning a period of 4 weeks.

Antiviral Therapy

We simulated therapy with a hypothetical antiviral agent that has the same mechanism of action as Remdesivir,[61] i.e. interference with RNA-dependent RNA polymerase. For this, 200 mg loading dose (≡initial plasma concentration C0 = 5 μM, assuming a molecular weight of 800 g/mol and a volume of distribution of 50 L), followed by 100 mg daily maintenance doses for 9 additional days via intravenous injection were simulated to mimic the pharmacological intervention. Treatment was started at one of the two time points previously mentioned. As a simplification, assuming one-compartment pharmacokinetics for the hypothetical drug with elimination rate constant kCl = 1 d–1, we use the plasma concentration kinetics as a surrogate for tissue concentration kinetics of the drug. The plasma concentration αV(t) is thus given bywhere .

Hence, the injection times are i = t, t + 1, ..., t + 9, and we define this set of times as S. The antiviral agent acts by inhibiting the production rate P of virus from infected cells in the body, such that P(t) is scaled by , where EC50 is the half maximal effective concentration of the drug, chosen to be 0.77 μM. Plasma concentration kinetics αV(t) of the drug are shown in the inset of Figure a for the scenario when the treatment was initiated on the day of onset of symptoms.

Figure 7

Simulated pharmacological interventions. Viral load kinetics in URT, LRT, and plasma following (a–c) treatment with a hypothetical antiviral agent, (d–f) treatment with interferon beta-1a, (g–i) and combination of the antiviral and interferon beta-1a, starting on the day of onset of symptoms (orange line) and 5 days post onset of symptoms (green), compared against no treatment (burgundy). Insets in panels a, d, and g show plasma concentration kinetics of the pharmacological agent/s for the treatment regimen that begins on the day of onset of symptoms. A lower bound at 102 copies mL–1 is imposed to represent the detectable limit of viral load in the body (black horizontal line); once the viral load goes below the detection limit, a vertical line is used to indicate the time of occurrence of this event.

As shown in Figure a–c, under the considered therapeutic regimen, the antiviral therapy shows a significant impact on viral load reduction compared to the no treatment scenario, irrespective of the timing of therapy initiation. However, early initiation of antiviral therapy led to a lesser total viral load and a faster reduction of the load, which is consistent with results of other studies.[19]

Interferon Therapy

We simulated treatment with interferon beta-1a, with and without the hypothetical antiviral agent (discussed above). For this, four subcutaneous injections (Dose = 44 μg each) of interferon beta-1a were administered every other day starting at one of the two time points previously mentioned. The plasma concentration kinetics of interferon beta-1a αI(t) is obtained by solving the following equation:where .

Hence, the injection times are i = t, t + 2, t + 4, t + 6, and we define this set of times as Q. The plasma concentration kinetics αI(t) of interferon beta-1a is shown in the inset of Figure d for the scenario when the treatment was initiated on the day of onset of symptoms. Here, the parameters kabsorp, Vsc, and kexcr, which represent the absorption of interferon beta-1a from the site of injection to the bloodstream, volume of the injection site, and excretion rate of interferon beta-1a, respectively, were estimated by fitting the above equation to literature derived concentration kinetics data.[62] The estimated parameter values are: kabsorp = 7.04 d–1,Vsc = 3.43 mL, and kexcr = 0.074 d–1. Of note, while interferon beta-1a acts in the same way as endogenous IFN, i.e. by inhibiting the production rate P of virus from infected cells, we model them separately to accommodate the possibility of unique pharmacokinetic properties of the two agents. Therefore, the first term of viral load kinetics in the ODE for organ i (see Methods) now becomes .

As shown in Figure d–f, interferon beta-1a therapy significantly reduces the viral load, and the effect is comparable to the hypothetical antiviral agent. Also, early initiation of therapy leads to a lesser total viral load and a faster reduction in viral load.

Combination Therapy

Finally, we tested the effect of the combination of the hypothetical antiviral agent and interferon beta-1a on viral load kinetics (plasma concentration kinetics in the inset of Figure g). As shown in Figure g–i, the combination therapy leads to a faster reduction in viral load compared to the effect of antiviral agent or interferon beta-1a alone, such that the viral load seems to fall below the detection threshold 2–3 days earlier. However, the need for clinical validation of this model prediction must be emphasized.

Conclusion

In summary, we have developed a semimechanistic (i.e., partially mechanistic and partially phenomenological) mathematical model that predicts the whole-body viral distribution kinetics of SARS-CoV-2 by incorporating cellular-scale viral dynamics, relevant physiological processes of viral transport, innate and adaptive immune response, and viral excretion. The model is well calibrated with published in vivo and clinical data and provides insights into the importance of various components of the immune system in controlling infection. Through global sensitivity analysis, we identified the key mechanisms that control infection and can be used as potential therapeutic targets for pharmacological intervention. Finally, we tested the potential of such therapeutic targets by simulating clinically relevant treatment options and identified the importance of the timing of treatment initiation and the effects of various therapies that may effectively suppress infection. As a limitation of the model, it is important to note that due to lack of clinical data for the viral load kinetics in extrapulmonary compartments, we had to rely on the transport parameters estimated from in vivo data, and the model predictions for the extrapulmonary compartments could not be clinically validated. Also, due to limited availability of data for cytokines in the hamster and clinical studies, the IFN-mediated inhibitory effects on viral production were calibrated with data for interferon-γ only. As clinical knowledge of the disease and its mechanisms improves, we will continue to fine-tune the model and integrate it with a physiologically based pharmacokinetic model to build an in-silico platform that can be used to simulate disease progression and a more complete pharmacokinetics of various test drugs and novel formulations.

Methods

Model Development

In each organ compartment, we consider the kinetics of three populations. The first is the concentration of healthy cells, denoted by y (t), where i is an arbitrary compartment. The second is the concentration of infected cells, denoted by y (t), and the third is the concentration of viral particles v(t). The compartments that we consider are i = [U, L, G, M, H, K, B], where U is the upper respiratory tract (URT), L is the lower respiratory tract (LRT), G is the gastrointestinal tract (GI), M is the mononuclear phagocytic system (MPS), which encompasses the liver and the spleen, H is the heart, K denotes the kidneys, and B is the brain. All organ compartments are connected via the plasma compartment P(t), in which we only consider the concentration of viral particles in systemic circulation.

We assume that, due to the time scale under consideration, the population of healthy cells does not replenish during the simulation window. Thus, the rate of change of healthy cells is modeled as a decreasing function, which depends on the interaction between healthy cells and the surrounding viral particles that infect them. The conversion rate of healthy cells to infected cells is denoted by the infection rate I.

The rate of change of total infected cell population within an organ is a function of three mechanisms. The first is an increase due to freshly infected healthy cells by viral particles, the second is death due to the cytopathic effects of viral infection, proportional to the number of infected cells and characterized by the rate constant D1, and the third is due to the interaction between effector CD8+ cells (concentration denoted by CD8*(t)) and the infected cells. The cytotoxic effect is measured in terms of the death rate constant DCD8. It is important to note that in the reduced form of the model, infected cell death due to effector CD8* cells is not included. In our model description, all activated immune cells are indicated with an asterisk (e.g., CD8*), while inactive or naïve populations are indicated by standard naming conventions (e.g., CD8+).

The concentration of viral particles in each organ compartment is influenced by different factors, some of which are organ-specific. However, in all organ compartments, the rate of change of viral particles is proportional to the number of infected cells with virus production constant P but is inversely proportional to the concentration of IFN, denoted IFN(t). IFN is part of the innate immune response that acts by suppressing the viral production rate of infected cells and is controlled by the effectiveness constant ε. In all compartments, viral particles may be neutralized due to the aggregation of antibodies on their surface proteins. The antibody concentration is represented by Ab(t), and the destruction rate of viral particles due to antibodies is denoted by DAb. Of note, neutralization of virus by antibodies is not incorporated in the reduced form of the model. Lastly, the viral load can also be reduced by tissue resident macrophages in the lungs (alveolar macrophages) and MPS (Kupffer cells and splenic macrophages) that engulf the viral particles and destroy them. The rate of removal of viral particles by macrophages is given by λMϕ.

We introduce the system of ordinary differential equations that governs the concentration kinetics as follows. Mechanisms particular to each compartment are discussed after each set of equations.

Equations for the URT:Note that the rate of change of viral particles in the URT depends on two clearance mechanisms. These consist of the migration of viral particles from the URT to the LRT and from the URT to the GI tract. The corresponding transport coefficients are C1 and C2, respectively.

Equations for the LRT:Viral particles in the LRT are received from the URT but are also lost to the systemic circulation, which is quantified through the pulmonary absorption coefficient A1.

Equation for IFN:As mentioned earlier, IFN limits the production of viral particles by infected cells. The rate of production of this cytokine is regulated by two effects. One is a zeroth order generation term GIFN, and the second is proportional to the cumulative population of infected cells in the respiratory tract. The proportionality constant in the second mechanism is PIFN. Further, IFN can degrade over time at a rate given by DIFN.

Equations for the GI tract:Within the GI tract, in addition to the standard mechanisms, viral particles can be lost to the liver via the hepatic portal vein. This rate is quantified through the intestinal absorption rate A2. Additionally, as discussed above, viral particles from the URT are transported to the GI tract at a rate C2.

Equations for the MPS:The rate of change of viral particles within the MPS can increase due to the particles collected from the GI tract (A2) and those incoming from systemic circulation through the hepatic artery (H). At the same time, particles can leave this compartment through the hepatic vein (H) or through the hepatobiliary excretion mechanism (E).

Equations for the heart:Similar to the MPS, the rate of change of viral particles in the heart depend on two fluxes. One is the incoming viral load from systemic circulation through the coronary arteries (CO), and the second is the outgoing viral particles via the coronary veins (CO).

Equations for the kidneys:Given that the size of the SARS-CoV-2 virus is ∼100 nm, we assume that renal excretion is not feasible.[63] Therefore, the viral kinetics in kidneys is dependent on the standard mechanisms of influx via the renal arteries (R) and outflux via the renal veins (R).

Equations for the brain:While the blood–brain barrier may be deterrent to the establishment of infection in the brain, we have included the brain compartment due to the lack of evidence for the former.[64] Analogous to the heart and kidney compartments, the brain receives viral particles delivered through systemic circulation via the internal carotid arteries (CA). Subsequently, viral particles can rejoin the plasma compartment by means of the internal jugular veins (J).

Equations for plasma:The equation for the plasma compartment incorporates all the outgoing fluxes of viral particles through arteries (H, R, CO, CA) and the incoming fluxes via veins (H, R, CO, j). Also, the viral load from LRT in eq gets added to plasma at rate A1. Lastly, similar to all the previous compartments, viral particles can be neutralized by antibodies at a rate DAb (in the full model only).

The model up to this point is the reduced model, and the only equation that captures the immune system is eq . This equation describes the change of concentration of IFN, which is part of the innate immune system. However, the adaptive immune system should activate to properly mount a full immune response. The equations that follow provide the remaining elements to initiate and maintain a humoral and cell-mediated adaptive immune response to make up the full model.

Equations –30 model the concentration of naïve antigen presenting cells (APCs) APC(t) in a given compartment i. These cells predominantly act as carriers, transporting remnants of viral particles to the lymphatic compartment to raise the alarm in the adaptive arm of the immune system. In order to keep the population of APCs at steady state, a replenishing mechanism is included. It is proportional to the difference between the current concentration APC(t), and the generic steady state value APC, divided by the absolute value of the same difference plus a modulating constant that we set = 1, i.e., . The rationale behind the previous definition is that if the concentration APC(t) is close to APC, the term is essentially zero, whereas if the difference is large, the denominator is equally large, and the ratio is close to one. This allows a near constant modulation rate that only acts whenever the population of macrophages is far from equilibrium. This reequilibration rate is denoted by GAPC. The APC population is also impacted by the interaction between APCs and invading viral particles. This causes the population of APCs to become activated at a rate TAPC, which is proportional to the product of the concentration of these cells and the viral load.

Equations for URT naïve APCs:Equations for the LRT naïve APCs:GI tract naïve APCs:MPS naïve APCs:Heart naïve APCs:Kidney naïve APCs:Brain naïve APCs:Lymphatic activated APCs:As mentioned earlier, the concentration of activated APCs, denoted by APC*(t), increases at a rate proportional to the conversion rate TAPC of naïve APCs to the activated state. At the same time, activated APCs can become incapacitated or simply eliminated after a certain period of time. These processes are pooled into the APC death rate constant DAPC.

The remaining equations in the model describe either lymphocyte or antibody concentrations. The lymphocytes under consideration are CD8+, CD4+, and B cells. The concentration of the naïve or inactivated version of these populations is represented by CD8(t), CD4(t), and BC(t), respectively. Similar to the APC case, the steady state concentration for each of these cell types is denoted with a horizontal bar on top of each variable. This value is maintained using a mechanism analogous to the APC case, i.e., through the term , where X is the corresponding cell type.

Furthermore, each of these lymphocytes can be promoted to its activated state, which we denote with an asterisk, X*. This conversion takes place when activated APCs interact with naïve lymphocytes, presenting them fragments of the ingested viral particles. The transformation rate from the naïve state to the activated one is represented by the constant T, where X indicates the cell type. For the case of CD4+ and CD8+, we use the common rate TT. The concentration of activated lymphocytes can fluctuate due to two mechanisms. One is the conversion of naïve cells, which was already discussed. The second is due to death caused by a variety of factors, such as exhaustion or apoptosis, and is expressed through the death constant D. B cells, however, follow a different mechanism, which we discuss later.

Naïve CD8:Active CD8*:Naïve B cells:Activated B cells:Activated B cells alone cannot neutralize the viral load. They must first transform into plasma cells, which then in turn produce antibodies that can carry out the viral neutralization. The transition from activated B cell to plasma cell is enabled by activated CD4* cells. This transformation can generate two types of plasma cells: short-lived (S) and long-lived (L). The transition rates are represented by the constants TBCS and TBCL, and the corresponding concentrations of plasma cells are denoted by PS(t) and PL(t), respectively. Lastly, each type possesses a particular death rate, indicated by the constant DPS for the short-lived type and TPL for the long-lived type.

Naïve CD4:Activated CD4:Short-lived plasma cells:Long-lived plasma cells:Antibodies:The final equation in the model describes the rate of change of antibody concentration, Ab(t). The sole contributors to the production of antibodies are short-lived and long-lived plasma cells. The corresponding production rates are given by GPS and GPL, respectively. We consider two mechanisms by which antibodies can be consumed or eliminated. One is the loss of antibodies due to their interaction with viral particles in a given compartment with rate constant λAb. The second is the elimination of antibodies by other factors independent of the viral load, e.g., degradation or clearance from tissues, which occurs at rate constant ClAb. A summary of the parameters used in the model is given in Tables –3.

Once the full set of parameters for the model was defined and fixed for the conditions corresponding to Figure , we proceeded to identify the relative effect of each parameter on the viral kinetics. This was done by systematically perturbing multiple parameters simultaneously and comparing the outcomes between the perturbed state and the original state. To quantify the differences between the two states, we introduced the following six criteria. The first three refer to the area under the curve (AUC) of the viral load curves in the URT, the LRT, and the plasma compartment for a period of 30 days. Mathematically,The next three denote the total time required for the viral load in the URT, LRT, and plasma compartment to reduce to <102 viral copies per mL, i.e., 2 in log base 10 units. In the next subsection, we compare each of the 6 criteria using global sensitivity analysis (GSA), applied to the 31 parameters.

Global Sensitivity Analysis

Local sensitivity analysis only allows one parameter to change at a time and thus may not shed light on the effect of parameter interactions in complex biological systems. Thus, we performed global sensitivity analysis (GSA;[39,65−67] i.e., multiple parameters can be varied simultaneously) to obtain a further understanding of the relationship between parameters and their combined effects on the model outputs of interest. Specifically, in our analysis, all parameters were perturbed at the same time, and we assumed a uniform distribution of values for each parameter. The range of values of each sample was PR ± 99% of PR, where PR is the reference value of a given parameter. To ensure comprehensive sampling of the multidimensional parameter space, we used a Latin hypercube sampling scheme. Subsequently, each sample (made up of 31 parameter values) was used to compute the kinetics, from which we extracted the 6 criteria mentioned above. In total, 10 batches of 5000 samples each were computed. Then, each batch was subjected to a multivariate linear regression analysis to obtain the regression coefficients of the parameters, and then averaging the regression coefficients of each parameter over 10 batches, we obtained the mean and variance of the sensitivity index (SI) for each parameter. Lastly, the sensitivity indices were ranked by means of a one-way ANOVA and a Tukey’s test. The higher the SI, the more significance a parameter holds. A detailed description of the GSA workflow can be found in refs (39 and 65−67). All analyses were performed in MATLAB R2018a.