Dynamically modeling the effective range of IL-2 dosage in the treatment of systemic lupus erythematosus.

Systemic lupus erythematosus (SLE) is a complex systemic autoimmune disease characterized by an overactive immune response to self-antigen. The overactivation of CD4+ Foxp3- conventional T cells (Tcons) and the inactivation of CD4+ CD25+ Foxp3+ regulatory T cells (Tregs) play important roles in the progression of SLE. Clinical trials showed that low-dose interleukin-2 (IL-2) is effective in treating SLE. Here, we developed a mathematical model involving Tcons, Tregs, natural killer (NK) cells, and IL-2 to simulate the dynamic processes involved in the treatment of SLE. We found an effective range of IL-2 dosage defined by the Tcon/Treg ratio in SLE treatment, termed the IL-2 dosage therapeutic window (IDTW). Our results showed that high levels of self-antigen result in a narrow IDTW and high post-treatment Tcon/Treg ratio. Furthermore, we proposed a classification method based on the ratio of pre-treatment Treg to CD4+ T cells to predict the treatment outcome of SLE patients.

Introduction

Systemic lupus erythematosus (SLE) is a non-organ-specific, complex, and chronic systemic autoimmune rheumatic disease (Lisnevskaia et al., 2014; Tsokos, 2011) characterized by the immune imbalance of T cell subsets, in particular, the downregulation of CD4+ CD25+ Foxp3+ regulatory T cells (Treg) and the upregulation of pathogenic CD4+ Foxp3− conventional T cells (Tcon), such as follicular helper T cells (Tfh) and T helper 17 cells (Th17) (Comte et al., 2016; He et al., 2013; Rose et al., 2019). This imbalance is caused by excessive immune response induced by persistent stimulation of self-antigen, such as double-stranded DNA (dsDNA) (Dean et al., 2000). The production of self-antigen in SLE patients is mainly attributed to defective apoptotic clearance (Shao and Cohen, 2011). Phagocytes fail to efficiently remove apoptotic material, which results in the capture of dsDNA fragments by antigen-presenting cells (APCs). In response to stimulation by self-antigen presented by APCs, Tcon is overactivated, suggesting the disease status of SLE patients. Previous studies showed that there is an imbalanced Treg/Tcon in SLE, especially Treg/Tfh. Treg/Tfh ratio was inversely correlated with SLE disease activity (SLEDAI) and proinflammatory cytokine levels (He et al., 2016; Miao et al., 2021), which suggested Tcon subsets contribute to the aggravation of SLE. In the clinical practice of SLE treatment, Tcon cells promote the inflammatory response, whereas Tregs have anti-inflammatory effects (Abbas et al., 2018). However, we still lack mechanistic understanding and clinical data to quantify the effects of Tcon and Treg on the inflammatory response. At this point, using the Tcon/Treg ratio to characterize SLE disease progression may be an appropriate approach at the current stage.

The treatment of SLE is challenging. Over the past decades, the treatment of SLE has relied on the administration of glucocorticoids and immunosuppressive agents (Tsokos, 2011). Recently, some studies reported that low-dose exogenous interleukin-2 (IL-2) supplementation can be a promising therapeutic strategy for several autoimmune-related diseases (Koreth et al., 2011; Saadoun et al., 2011), including SLE (He et al., 2016; Spee-Mayer et al., 2016). IL-2 is essential to maintain the functionality of Tregs and restore immune homeostasis (Dooms and Abbas, 2010; Humrich et al., 2010). Most applications of low-dose IL-2 in SLE treatment have been reported as effective and well-tolerated. However, high-dose IL-2 can lead to adverse events, and IL-2-induced natural killer (NK) cell proliferation has been reported to correlate with the incidence and severity of adverse events (Humrich et al., 2019). Owing to the limitations of clinical use, the role and mechanism of high-dose IL-2 in SLE patients have still not been fully elucidated (Huang and Stott, 1995). Further phase I/II clinical trials are underway to provide more detail on the efficacy and safety of IL-2 therapy (He et al., 2020; Klatzmann and Abbas, 2015).

Here, we established a mathematical model to reveal the dynamics of Tcon, Treg, and NK populations during IL-2 treatment in SLE patients. We found an effective IL-2 range in treating SLE patients defined by the Tcon/Treg ratio, termed the IL-2 dosage therapeutic window (IDTW), and we analyzed the formation processes of IDTW. Furthermore, we found that the width of IDTW and treatment effect depends on self-antigen levels that usually vary widely among different patients. Considering the heterogeneity of SLE patients, we proposed a classification method based on pre-treatment Treg in CD4+ T cell ratio to divide SLE patients into IL-2-responsive and IL-2-unresponsive two groups before treatment, which may be used to assess the ultimate benefit of treating SLE patients with IL-2. Finally, we built a two-variable toy model and analyzed the nullclines in phase space to reveal the underlying mechanism leading to treatment failure with high-dose IL-2. Our mathematical models provide deeper insights into the efficacy of IL-2 therapy and may further improve current treatment strategies.

Results

Network in SLE treatment and the mathematical model

The background of IL-2 therapy in SLE patients

As SLE is a complex non-organ-specific disease and our focus is on IL-2 treatment, we only considered the immune response in peripheral blood. During the SLE disease process, self-antigens stimulate APCs to proliferate and present specific antigens to naive T (TN) cells, which triggers TNs to differentiate into Tcons and Tregs. Here, Tcon cells contain mature CD4+ cells other than Tregs, including Th1, Th2, Th17, Tfh, etc. Activated Tregs suppress the overactivation of Tcons and the function of APCs to avoid immune imbalance (Grover et al., 2021). However, in SLE patients, Tregs are shown to be functionally deficient, thus losing their capacity to regulate Tcon activity (Vitales-Noyola et al., 2017). Previous research has proved that the lack of IL-2 secretion capacity in SLE patients is associated with the defects of function in Tregs (Boyman et al., 2015; Karampetsou et al., 2016; Whitehouse et al., 2017). The administration of exogenous low-dose IL-2 could restore the functionality of Treg (Spee-Mayer et al., 2016).

Tcon and TN cells can secrete IL-2 and release IL-2 extracellular. In Figure 1C, we quantitatively assessed the endogenous IL-2 secretion capacity in SLE patients, we tested the proportion of IL-2 producing CD4+ effector T cells (Teffs), Tregs, CD8+ T cells, double-negative T cells (DNT), NK cells, and natural killer T cells (NKT) in peripheral blood of eight SLE patients and five healthy controls (Figure 1C). CD4+ Teffs are composed of Tcon and naive CD4+ T cells without Tregs (Figure S3), and NK cells correspond to NK in Figure 1A. We found that the ratio of IL-2 producing CD4+ T cells in SLE patients was less than that in healthy controls, suggesting the weak capacity of SLE patients for endogenous IL-2 production. Thus, we substantially downregulated the IL-2 secretion capacity of TNs and Tcons in the SLE model compared with healthy controls (we compared the healthy person with SLE patients in Figure S8). No difference was observed in the IL-2 expression by NKT cells and DNT cells between SLE patients and healthy controls. In addition, when administering IL-2 to treat SLE patients, exogenous IL-2 is present in peripheral blood.

Figure 1

Schematic diagrams of IL-2 related interaction network in SLE patient (A) and IL-2 releasing in peripheral blood (B)

(A) Upon Ag stimulation presented by antigen-presenting cells (APCs), naive T (TN) cells differentiate into Tcons and Tregs. Tcons and TNs secrete IL-2 that can be captured by IL-2 receptors on Tcons and Tregs. SLE patients show a functional deficit of Tregs and weak endogenous IL-2 secretion capacity (dashed lines). Administration of exogenous IL-2 can restore Treg functionality to suppress the overactivation of Tcons.

(B) In peripheral blood with a total volume of , the free volume in peripheral blood equals minus the volume of all cells. The release of IL-2 by TNs and Tcons involves a volumetric change from intracellular to extracellular.

(C) Percentages of IL-2-producing CD4+ Teff, double-negative T (DNT) cells, natural killer (NK) cells, and natural killer T (NKT) cells in eight SLE patients and five healthy controls were obtained by flow cytometry. Horizontal lines in scatterplots indicate the median and interquartile range, respectively. The Mann–Whitney U test was performed for non-parametric data.

Tcons, Tregs, and NK cells compete for IL-2 through the IL-2 receptor (IL-2R) on their surface, IL-2R binds IL-2 to form the IL-2/IL-2R complex and thereby mediates the proliferation of these cells. As the density of IL-2R on Tregs is higher than on Tcons (Abbas et al., 2018), IL-2 is preferentially captured by the IL-2R constitutively expressed on Treg (ichi Matsuoka et al., 2013) and forms the IL-2/IL-2R complex. Tcon, despite the low density of surface IL-2R, also competes for IL-2 to form the IL-2/IL-2R complex, constituting a positive feedback (Busse et al., 2010). NK cells, approximately 90% are CD56 NK cells (He et al., 2020; Spee-Mayer et al., 2016), mostly express the dimeric IL-2R (Abbas et al., 2018), which also competes for IL-2 and forms the IL-2/IL-2R. Besides, CD25 is the chain of IL-2R, owing to soluble CD25 (sCD25) levels that were significantly higher in SLE patients compared with the healthy controls (El Shafey et al., 2008), we also considered the binding of sCD25 to IL-2 in our SLE model.

The network and mathematical model of SLE treatment

To reveal the mechanism of IL-2 therapy in SLE patients, in Figure 1A, we proposed a simplified network involving Tcons, Tregs, TN cells, APCs, NK cells, and IL-2, where IL-2 is secreted by Tcons and TNs and competed by IL-2R on the surface of Tcons, Tregs, and NK cells. As Tregs are shown to be functionally deficient, thus losing their capacity to regulate Tcon activity (Vitales-Noyola et al., 2017), we assumed in our model that Tregs in SLE patients can be divided into functional Treg (Treg) and dysfunctional Treg (Treg) to simulate the functional defects in Tregs accordingly. The administration of exogenous low-dose IL-2 could transform Treg to Treg and restore the functionality of Treg (Spee-Mayer et al., 2016).

In Figure 1B, we illustrated the IL-2 releasing process. We set the total volume in peripheral blood as , while the volumes of all APCs, TNs, NKs, Tcons and Tregs are , , , , and , respectively. After IL-2 is released from the cells, it is uniformly distributed in the volume of , where denotes the volume of free space in peripheral blood and is equal to minus the total volume of all cells in .

Based on the analysis, we developed an ordinary differential equations (ODEs) model to investigate the dynamic processes of the immune response during IL-2 treatment in our model (Khailaie et al., 2020; Sontag, 2017). , , , , and denote the density of APC, TN, Treg, Treg, and NK cell, respectively, they all are in the unit mm−3. denotes the concentration of extracellular IL-2 with the unit nM. , , and represent the density of IL-2/IL-2R complex on the surface of Tcon, Treg, and NK cell respectively, all they are in the unit m−2. The dynamic equations are shown as follows:

The meaning and value of each parameter in the equations can be found in Table 1, and the initial value of all variables can be found in Table 2. We described the equations in three parts: cell proliferation (Equations 1, 2, 3, 4, 5, and 6), extracellular IL-2 (Equation 7), and IL-2 binding with IL-2R (Equations 8, 9, and 10) as follows.where denotes that APCs have an environmental carrying capacity . denotes the proliferation rate of APC and is the proliferation coefficient of APCs induced by antigen. We used hill equation to describe the cell proliferation and chose the order in our model owing to the complex interactions between variables (different choices of can be found in Figure S6). The last term describes the suppression of APCs by Treg with a suppression rate .

Table 1

The descriptions, values, and sources of kinetic parameters in our model

Parameters	Descriptions	Values	Sources
p_APC	Renewal rate of APCs	1 mm3h−1	Est. from (Friedl and Gunzer, 2001)
p_TN	Renewal rate of TN cells	20 mm3h−1	Est. from (Bains et al., 2009)
p_Tcon	Maximum differentiation rate from TNs to Tcons induced by APCs	0.018 h−1	Est. from (Bains et al., 2009)
p_Tregd	Maximum differentiation rate from TNs to Tregs induced by APCs	0.002 h−1	(He et al., 2016)
p_NK	Renewal rate of NK cells	0.6 mm3h−1	(Humrich et al., 2019)
k_APC	Proliferation rate of APCs	0.05 h−1	(Burroughs et al., 2006)
k_Tcon	Proliferation rate of Tcons	0.05 h−1	(Burroughs et al., 2006)
k_Tregd	Proliferation rate of Tregd	0.05 h−1	(Burroughs et al., 2006)
k_Tregf	Proliferation rate of Tregf	0.05 h−1	(Burroughs et al., 2006)
k_NK	Proliferation rate of NK cells	0.05 h−1	(Burroughs et al., 2006)
k_trans	Transformation rate from Tregd to Tregf	0.05 h−1	Est.
d_APC	Apoptosis rate of APCs	0.01 h−1	(Price et al., 2015)
d_TN	Apoptosis rate of TNs	0.01 h−1	(Price et al., 2015)
d_Tcon	Apoptosis rate of Tcons	0.01 h−1	(Price et al., 2015)
d_Tregd	Apoptosis rate of Tregd	0.01 h−1	(Price et al., 2015)
d_Tregf	Apoptosis rate of Tregf	0.01 h−1	(Price et al., 2015)
d_NK	Apoptosis rate of NK cells	0.01 h−1	(Price et al., 2015)
d_trans	Transformation rate from Tregf to Tregd	0.01 h−1	Est.
K_APC	Environmental carrying capacity of APCs	200 mm−3	Est. from (Friedl and Gunzer, 2001)
K_Tcon	Environmental carrying capacity of Tcons	3000 mm−3	(Shete et al., 2010)
K_Treg	Environmental carrying capacity of all Tregs	400 mm−3	(He et al., 2016)
K_NK	Environmental carrying capacity of NK cells	400 mm−3	(Humrich et al., 2019)
S_Tregf,APC	Suppression rate of Tregf to APCs	0.03 h−1	Est.
S_Tregf,Tcon	Suppression rate of Tregf to Tcons	0.02 h−1	Est.
Ag	Concentration of self-antigen	5 A.U.	Est.
K_Ag,APC	Proliferation coefficient of APCs induced by antigens	1 A.U.	Est.
K_APC,TN	Differentiation coefficient of TNs induced by APCs	100 mm3	Est.
K_Tregf,APC	Suppression coefficient of APCs induced by Tregf	50 mm3	Est.
K_Tregf,Tcon	Suppression coefficient of Tcons induced by Tregf	500 mm3	Est.
K_trans	Activation coefficient of Tregf induced by IL-2/IL-2R complex	15 μm2	Est.
R_Tcon	Proliferation coefficient of Tcons induced by IL-2/IL-2R complex	0.25 μm2	Est.
R_Tregd	Proliferation coefficient of Tregd induced by IL-2/IL-2R complex	15 μm2	Est.
R_Tregf	Proliferation coefficient of Tregf induced by IL-2/IL-2R complex	15 μm2	Est.
R_NK	Proliferation coefficient of NK cells induced by IL-2/IL-2R complex	0.25 μm2	Est.
n	The order of hill equations	2
k_Tcon,IL-2	Secretion rate of IL-2 by Tcons	1 nMh−1	Est. from (Busse et al., 2010)
k_TN,IL-2	Secretion rate of IL-2 by TNs	1 nMh−1	Est. from (Busse et al., 2010)
ε	Conversion factor	1 μm2nM	Est.
δ	Volume effect factor	5 × 10−7 mm3	Est.
k_αβγ	Binding rate of IL-2 and IL-2Rαβγ on Tcons and Tregs	111.6 nM−1h−1	(Busse et al., 2010)
k_α	Binding rate of IL-2 and IL-2Rα	k_αβγ/3 nM−1h−1	Est. from (Busse et al., 2010)
k_βγ	Binding rate of IL-2 and IL-2Rβγ on NK	k_αβγ/3 nM−1h−1	Est. from (Busse et al., 2010)
c_αβγ	Dissociation rate of IL-2/IL-2Rαβγ complex	0.83 h−1	(Busse et al., 2010)
c_βγ	Dissociation rate of IL-2/IL-2Rβγ complex	0.83 h−1	(Busse et al., 2010)
d_αβγ	Degradation rate of IL-2/IL-2Rαβγ complex	1.7 h−1	(Busse et al., 2010)
d_βγ	Degradation rate of IL-2/IL-2Rβγ complex	1.7 h−1	(Busse et al., 2010)
d_IL-2	Degradation rate of extracellular IL-2	0.1 h−1	(Busse et al., 2010)
IL-2_exo	Concentration of exogenous IL-2	0.02 nMh−1	Est. from (He et al., 2016)
IL-2RNKtot	Average surface density of all IL-2R on NK cells	0.3 μm−2	Est. from (Busse et al., 2010)
IL-2RTcontot	Average surface density of all IL-2R on Tcons	0.3 μm−2	Est. from (Busse et al., 2010)
IL-2RTregtot	Average surface density of all IL-2R on Tregs	30 μm−2	Est. from (Busse et al., 2010)
[sCD25]	Average concentration of sCD25	0.01 nM	Est. from (El Shafey et al., 2008)

Est. is the abbreviation of estimated.

Table 2

The initial value of variables in our model, corresponding to the pre-treatment state of SLE patient

Variables	Values	Units	Variables	Values	Units
[APC]	178	mm−3	[TN]	793	mm−3
[Tcon]	1098	mm−3	[Tregd]	122	mm−3
[Tregf]	5	mm−3	[NK]	60	mm−3
[IL-2]	0.0015	nM	[IL-2/IL-2R_Tcon]	0.019	μm−2
[IL-2/IL-2R_Treg]	1.9	μm−2	[IL-2/IL-2R_NK]	0.0066	μm−2

This part motivates the equations describing the dynamic of cell density. Equation 1 describes the dynamic of APC cell density, where denotes the renewal rate and denotes the apoptosis rate of APCs. Under antigen stimulation , we described the proliferation of APC cells as:

Equation 2 describes the dynamic of naive T cell density, where denotes the renewal rate and denotes the apoptosis rate of TNs. TN cells will differentiate into Tcon and Treg cells with the differentiation rate and under stimulation by APC cells.

Equation 3 describes the dynamic of Tcon cell density, where denotes the renewal rate and denotes the apoptosis rate of Tcons. The average surface density of IL-2/IL-2R complex on Tcons is , and the number of IL-2/IL-2R complex on a single cell is ( is the radius of a single cell). We gave the proliferation term of Tcons induced by the IL-2/IL-2R complex with a proliferation rate aswhere is the environmental carrying capacity of Tcons. Dividing the equation above by , the proliferation term of Tcons can be simplified aswhere is the proliferation coefficient of Tcons induced by IL-2/IL-2R. The formation of IL-2/IL-2R will be discussed in the next chapter. Similar to APC cells, Tcons will also be suppressed by Treg with a suppression rate .

Equation 4 describes the dynamic of Treg cell density, where denotes the renewal rate and denotes the apoptosis rate. We described the proliferation of Treg induced by IL-2/IL-2R with a term similar to Tcons. Besides, the IL-2/IL-2R complex also induces the activation from Treg to Treg with an activation rate as

Equation 5 describes the dynamic of Treg cell density, which is activated from Treg cells induced by IL-2/IL-2R. The third term shows the proliferation of Treg induced by IL-2/IL-2R with a term similar to Tcons. The deactivation of Treg to Treg is given by .

Equation 6 describes the dynamic of NK cell density with a renewal rate of and an apoptosis rate of . The proliferation term of NK cells induced by IL-2/IL-2R is similar to those of Tcons.

In this part, we showed the dynamic of extracellular IL-2 concentration in Equation 7. As our model is cross-scale, when considering the process of IL-2 secretion or IL-2 binding, we introduced the volume effect factor to describe the secretion of IL-2 from intracellular to extracellular (see STAR Methods). The term describes IL-2 secreted by Tcons can be written as , where is the IL-2 secretion rate by Tcons. The term describes IL-2 secreted by TNs can be written as , where is the IL-2 secretion rate by TNs.

Then, we derived the increment of extracellular IL-2 by dissociation of the IL-2/IL-2R complex. Both binding of IL-2 to the IL-2R and dissociation of IL-2/IL-2R occur on the cell surface, we assumed that IL-2R and IL-2/IL-2R complex are evenly distributed intracellular by introducing conversion factor (see STAR Methods). Considering that Tcon cells express the trimer IL-2R , we set the dissociation rate of IL-2/IL-2R on Tcons as and the decrement of IL-2 concentration in Tcon in d is . Thus, extracellular IL-2 concentration increased by the dissociation from IL-2/IL-2R on all Tcon cells can be written as . Note that the binding of IL-2 with IL-2R on Tcons has a similar process to the dissociation process, we can write the extracellular IL-2 concentration decreased by the binding to IL-2R on all Tcons as , where is the binding rate of IL-2 and IL-2R on Tcons. Similarly, we can derive the changes of extracellular IL-2 concentration by binding to IL-2R or dissociation of IL-2/IL-2R on Tregs in the same way. As most NK cells express the dimeric IL-2R , we used , and to distinguish the binding, dissociation and degradation of IL-2R on NK cells from IL-2R on T cells. Soluble CD25 (IL-2R ) also bind to IL-2 as described by the term , where is the concentration of sCD25 and is the binding rate of IL-2 with IL-2R . As the mechanism of sCD25 production is not fully understood and the interaction between sCD25 and immune cell subpopulations is complex, we set the concentration of sCD25 as a parameter through the rapid equilibrium assumption (see STAR Methods). Last, represents exogenous IL-2 concentration during IL-2 therapy.

In the third part, we derived the formation of the IL-2/IL-2R complex. Equation 8 describes the dynamic of the average density of the IL-2/IL-2R complex on Tcons. IL-2R on Tcons binds to extracellular IL-2 to form the IL-2/IL-2R complex, denoted by in the model (Molecular perspective of IL-2 binding process can be found in STAR Methods). Besides, IL-2/IL-2R complex on Tcons will decrease by dissociation and degradation, as indicated by the terms and . In our model, we assumed that the sum of and is a constant .

Equations 9 and 10 denote the dynamic of the average density of IL-2/IL-2R complex on Treg and NK cells, which can be derived in the same way as IL-2/IL-2R on Tcons. Table 3 briefly shows the biological meanings and parameters of each term in the equations.

Table 3

The scheme shows all model components with involved parameters, biological meanings and derivations

Equation 1[APC]	1: p_APC	2: −dAPC[APC]	3: k_APCAg^nAg^n+KAg,APCn[APC](1−[APC]K_APC)	4: −S_Tregf,APC[APC]^n[APC]^n+KTregf,APCn[Tregf]
Parameters	p_APC	dAPC	k_APC, Ag, K_Ag,APC, K_APC, n	S_Tregf,APC, K_Tregf,APC, n
Meaning	Renew	Apoptosis	Proliferation	Suppression
Derivation	(1)	(1)	(1)	(1)
Equation 2[TN]	1: p_TN	2: −d_TN[TN]	3: −p_Tcon[APC]^n[APC]^n+KAPC,TNn[TN]	4: −p_Tregd[APC]^n[APC]^n+KAPC,TNn[TN]
Parameters	p_TN	d_TN	p_Tcon, K_APC,TN, n	p_Tregd, K_APC,TN, n
Meaning	Renew	Apoptosis	Differentiation	Differentiation
Derivation	(1)	(1)	(1)	(1)
Equation 3[Tcon]	1: p_Tcon[APC]^n[APC]^n+KAPC,TNn[TN]	2: −d_Tcon[Tcon]	3: k_Tcon[IL-2/IL-2R_Tcon]^n[IL-2/IL-2R_Tcon]^n+RTconn[Tcon](1−[Tcon]K_Tcon)	4: −S_Tregf,Tcon[APC]^n[APC]^n+KTregf,Tconn[Tregf]
Parameters	p_Tcon, K_APC,TN, n	d_Tcon	k_Tcon, R_Tcon, K_Tcon, n	S_Tregf,Tcon, K_Tregf,Tcon, n
Meaning	Differentiation	Apoptosis	Proliferation	Suppression
Derivation	(1)	(1)	(1)	(1)
Equation 4[Tregd]	1: p_Tregd[APC]^n[APC]^n+KAPC,TNn[TN]	2: −d_Tregd[Tcon]	3: k_Tregd[IL-2/IL-2R_Treg]^n[IL-2/IL-2R_Treg]^n+RTregdn[Tregd](1−[Treg]K_Treg)	4: −k_trans[IL-2/IL-2R_Treg]^n[IL-2/IL-2R_Treg]^n+Ktransn[Tregd]
Parameters	p_Tregd, K_APC,TN, n	d_Tregd	k_Tregd, R_Tregd, K_Treg, n	k_trans, K_trans, n
Meaning	Differentiation	Apoptosis	Proliferation	Transformation
Derivation	(1)	(1)	(1)	(1)
	5: d_Ttrans[Tregf]
Parameters	d_Ttrans
Meaning	Transformation
Derivation	(1)
Equation 5[Tregf]	1: k_trans[IL-2/IL-2R_Treg]^n[IL-2/IL-2R_Treg]^n+Ktransn[Tregd]	2: −d_Tregf[Tcon]	3: k_Tregf[IL-2/IL-2R_Treg]^n[IL-2/IL-2R_Treg]^n+RTregfn[Tregf](1−[Treg]K_Treg)	4: −d_Ttrans[Tregf]
Parameters	k_trans, K_trans, n	d_Tregf	k_Tregf, R_Tregf, K_Treg, n	d_Ttrans
Meaning	Transformation	Apoptosis	Proliferation	Transformation
Derivation	(1)	(1)	(1)	(1)
Equation 6[NK]	1: p_NK	2: −d_NK[NK]	3: k_NK[IL-2/IL-2R_NK]^n[IL-2/IL-2R_NK]^n+RNKn[NK](1−[NK]K_NK)
Parameters	p_NK	d_NK	k_NK, R_NK, K_NK, n
Meaning	Renew	Apoptosis	Proliferation
Derivation	(1)	(1)	(1)
Equation 7[IL-2]	1: k_Tcon,IL-2[Tcon]δ	2: k_TN,IL-2[TN]δ	3: −d_IL-2[IL-2]	4: −k_αβγεδ[IL-2][IL-2R_Tcon][Tcon]
Parameters	k_Tcon,IL-2, δ	k_TN,IL-2, δ	d_IL-2	k_αβγ, ε, δ
Meaning	Secretion	Secretion	Degradation	Binding
Derivation	(2)	(2)	(2)	(2)
	5: −k_αβγεδ[IL-2][IL-2R_Treg][Treg]	6: −k_βγεδ[IL-2][IL-2R_NK][NK]	7: c_αβγεδ[IL-2/IL-2R_Tcon][Tcon]	8: c_αβγεδ[IL-2/IL-2R_Treg][Treg]
Parameters	k_αβγ, ε, δ	k_βγ, ε, δ	c_αβγ, ε, δ	c_αβγ, ε, δ
Meaning	Binding	Binding	Dissociation	Dissociation
Derivation	(2)	(2)	(2)	(2)
	9: c_βγεδ[IL-2/IL-2R_NK][NK]	10: −k_α[IL-2][sCD25]	11: IL-2_exo
Parameters	c_βγ, ε, δ	k_α, [sCD25]	L-2_exo
Meaning	Dissociation	Binding	Injection
Derivation	(2)	(2)	(2)
Equation 8[IL-2/IL-2R_Tcon]	1: k_αβγ[IL-2][IL-2R_Tcon]	2: −c_αβγ[IL-2/IL-2R_Tcon]	3: −d_αβγ[IL-2/IL-2R_Tcon]
Parameters	k_αβγ	c_αβγ	d_αβγ
Meaning	Binding	Dissociation	Degradation
Derivation	(3)	(3)	(3)
Equation 9[IL-2/IL-2R_Treg]	1: k_αβγ[IL-2][IL-2R_Treg]	2: −c_αβγ[IL-2/IL-2R_Treg]	3: −d_αβγ[IL-2/IL-2R_Treg]
Parameters	k_αβγ	c_αβγ	d_αβγ
Meaning	Binding	Dissociation	Degradation
Derivation	(3)	(3)	(3)
Equation 10[IL-2/IL-2R_NK]	1: k_βγ[IL-2][IL-2R_NK]	2: −c_βγ[IL-2/IL-2R_NK]	3: −d_βγ[IL-2/IL-2R_NK]
Parameters	k_βγ	c_βγ	d_βγ
Meaning	Binding	Dissociation	Degradation
Derivation	(3)	(3)	(3)

(1), (2) and (3) denote the corresponding parts in Section “The network and mathematical model of SLE treatment”.

Our simulation results suggest an effective IL-2 range in SLE treatment

Through clinical efforts, administration of low-dose IL-2 has been proved effective in regulating the overactivation of Tcon in SLE patients (He et al., 2016; Spee-Mayer et al., 2016). However, owing to the limitations of clinical trials, it remains difficult to fully understand the role of IL-2 in SLE treatment, especially the effect of high-dose IL-2. We can analyze the immune response with different doses of exogenous IL-2 through our model.

Under adequate and persistent self-antigen stimulation, SLE patients usually have high Tcon and low Treg in CD4+ T cell ratios before treatment, representing the disease state. We simulated and illustrated the influence of four doses of IL-2 administration (0.01, 0.02, 0.05, 0.10 nM/h) on the density of Tcon, NK, and Treg cells (Figures 2A–2C). In the simulation, administration of 0.01 and 0.02 nM/h IL-2 have no effect on Tcon density, but slightly increases NK cell density and significantly increases Treg density, showing that the immune system has been effectively restored to homeostasis. The simulated changes in cell subpopulations caused by 0.01 and 0.02 nM/h IL-2 are qualitatively consistent with clinical experiments (He et al., 2016, 2020). However, the administration of 0.05 and 0.10 nM/h IL-2 doses can also significantly increase Tcon and NK cell density rather than only inducing Treg proliferation. The induction of Tcon and NK cell proliferation is often associated with higher severity of adverse events (Humrich et al., 2019). Furthermore, we showed the steady state of post-treatment Tcon, NK cell and Treg density in response to continuous IL-2 doses (Figures 2D–2F). These results suggest an effective IL-2 range for IL-2 therapy in SLE patients, as the density of IL-2R on the surface of Tregs is higher than on Tcons and NK cells, during this IL-2 range Tregs should activate and proliferate while NK and Tcon cells will not proliferate. In specific, we termed this range as the IDTW and selected the Tcon/Treg ratio to define the IDTW in our model. We demonstrated the post-treatment Tcon/Treg ratio in response to different IL-2 doses in Figure 2G. The pre-treatment Tcon/Treg ratio is 8.6, the post-treatment Tcon/Treg ratio is 3.5 under 0.01 nM/h IL-2 and the post-treatment Tcon/Treg ratio stabilize at 5.8 under 0.10 nM/h exogenous IL-2. Here, we defined the IDTW as reducing the Tcon/Treg ratio below 5.5 (see IDTW defined by other Tcon/Treg ratios in Figure S7; see parameters sensitivity analysis of IDTW in Figure S4). As shown in Figure 2G, the IDTW range is approximately between 0.004 and 0.045 nM/h, which approximately corresponds to 0.1 million and 1.125 million IU every day in clinical practice (see STAR Methods for the unit conversion). The IDTW range calculated by our model is close to the maximum IL-2 dose (1.5 million IU per day) estimated by Humrich et al. from clinical trials (Humrich et al., 2019).

Figure 2

Simulating the effective IL-2 range in treating SLE patients

The parameters in our simulation can be found in Table 1, and the initial value of variables can be found in Table 2.

(A–C) Time evolution of Tcon (A), NK cell (B), and Treg (C) with four exogenous IL-2 doses (0, 0.01, 0.02, 0.05, and 0.10 nM/h). With continuously administered IL-2, the system reaches a steady state after around 300 h.

(D–F) The post-treatment Tcon (D), NK cell (E), and Treg (F, pink line denotes the post-treatment Treg) to different IL-2 doses.

(G) The post-treatment Tcon/Treg ratio in response to different IL-2 doses (Ag = 5). Here the IDTW is defined to reduce Tcon/Treg below 5.5.

(H) For SLE patients with different antigen levels, we demonstrated the steady state of post-treatment Tcon/Treg ratio with various IL-2 doses (3D view in the left panel and 2D view in the right panel). The result shows that antigen level affects the range of IDTW and post-treatment Tcon/Treg ratio. The interval between the black line denotes the IDTW, which changes with patients’ self-antigen levels, and the white lines show the contour of the Tcon/Treg ratio. The value of the Tcon/Treg ratio is represented by the color bar.

We analyzed the mechanism of the formation of IDTW during IL-2 treatment. Owing to the constitutive and high expression of IL-2R on Treg, IL-2 in the IDTW is preferentially captured by Tregs and rescues the functionality of Treg. Tregs, in the functional state, effectively suppresses Tcons and restores homeostasis to the immune system. However, with an IL-2 dose higher than the IDTW, the capture, activation ,and then proliferation of Tcons is more than can be suppressed by Treg, and Tcon continues to proliferate. Thus, we can say that the overall formation of IDTW can be attributed to the greater competitiveness of Treg to IL-2 dose in the IDTW than Tcon through the constitutive expression of IL-2R.

Next, we investigated the effect of self-antigen level on the treatment process, especially the properties of IDTW, based on our mathematical model. In Figure 2H, we simulated the range of IDTW and the post-treatment Tcon/Treg ratio in SLE patients with different self-antigen levels. We found that the width of IDTW is inversely related to the self-antigen level of SLE patients, i.e., SLE patients with low self-antigen levels (Ag 1 in Figure 2H) have a wide range of IDTW. Besides, SLE patients with low self-antigen levels show a low post-treatment Tcon/Treg ratio with IL-2, demonstrating a good therapeutic effect. In comparison, the post-treatment Tcon/Treg ratio is relatively high in patients with high self-antigen levels when treated with IL-2 (Figure 3H). Considering the effects on IDTW width and the post-treatment Tcon/Treg ratio, we think that patients with lower self-antigen levels tend to have a better prognosis in IL-2 therapy compared with patients with high self-antigen levels.

Figure 3

Comparison of simulation results and clinical data

(A) We compared the evolution of the Treg ratio in CD4+ T cells under periodic administration of IL-2 in simulation with that given in He et al.’s clinical treatment data (He et al., 2016) (n = 23 patients, data are median with IQR). The red line denotes our simulated Treg ratio and the pink line denotes our simulated Treg ratio. The shaded parts represent the administration of IL-2 in SLE patients. They administered IL-2 (1 million IU, 0.02 nM/h) subcutaneously every other day for two weeks, and then followed by a two-week break. To compare time in the model with real time, we set  = .

(B) We compared the simulated evolution of the Tcon/Treg ratio with He et al.’s clinical treatment data (He et al., 2016) (n = 23 patients, data are median with IQR). The black solid line denotes our simulated Tcon/Treg ratio, and the black dotted line denote Tcon/Treg = 5.5. The blue line denotes our simulated Tcon ratio, and the value is shown on the right Y axis. The IL-2 therapy regime is the same as (A).

(C) We compared the NK cell ratio under periodic administration of IL-2 with He et al.’s data (He et al., 2020) (data are mean with SD). The meaning of symbols and the IL-2 administration regime is the same as (A). We set the number of lymphocytes to 890 mm−3 in this data set.

(D) shows the simulated pre- (pink line) and post-treatment (0.02 nM/h IL-2, red line) Tcon/Treg ratio, varying according to patients’ self-antigen level.

(E) IL-2 therapy reduces the Tcon/Treg ratios in SLE patients (n = 23 patients, p < 0.001 by paired t-test), and the bars denote the median of the Tcon/Treg ratios with IQR. The IL-2 therapy regime is the same as (A), and the IL-2 dose corresponds to 0.02 nM/h.

(F) We analyzed the correlation between dsDNA titers and post-treatment Tcon/Treg ratios in SLE patients (n = 23 patients). The Pearson correlation coefficient r is 0.58, and the p-value is 0.004 (the p-value for Pearson’s correlation coefficient uses the t-distribution).

Both computational results and clinical data suggest a decrease in the Tcon/Treg ratio after IL-2 treatment

Here, we semi-quantitatively compared the IL-2 trial results with simulation results through our model. Considering the interferences in clinical IL-2 trials, such as the combination of IL-2 with other drugs, the wide range of patient ages, and the limited sample size, the simulation results of our model will have some deviations from the clinical data.

He et al. treated SLE patients with low-dose recombinant human interleukin-2 (1 million IU rhIL-2) every other day for two weeks and followed by a two-week break off the drug (He et al., 2016). We introduced the data and compared the percentage of Treg cells in CD4+ T cells with our simulation result. By comparison, both the computational result and clinical trial demonstrate a significant increase of Treg in CD4+ T cell ratio during the IL-2 administration cycle and demonstrate a rapid decrease after discontinuation. We also simulated the change in Treg during the IL-2 administration cycle, which shows a similar trend to the Treg ratio. Besides, we calculated the changes in the Tcon/Treg ratio during IL-2 treatment and compared it with clinical data (He et al., 2016) in Figure 3B, both the results showed that IL-2 treatment successfully reduces Tcon/Treg below 5.5 and was therapeutically effective in SLE patients. Besides, He et al. tested the ratio of NK cells to lymphocytes and found an increase in NK cell ratio after IL-2 therapy (He et al., 2020). The regime and doses of IL-2 administration are the same as Figure 3A. We compared the simulation results of the NK cell ratio with clinical data in Figure 3C. Both the model and the experiment suggest that IL-2 therapy can lead to the proliferation of NK cells. However, the NK cell will rapidly decrease to normal levels after discontinuation of IL-2 (Figure 3C). We also showed the comparison of NK cell number density from Humrich et al. with our simulation (Humrich et al., 2019) in Figure S2, suggesting a similar trend with the NK ratio.

We examined the influence of antigen levels on treatment efficacy. Our simulation suggested two main results that 0.02 nM/h IL-2 can significantly reduce the Tcon/Treg ratio and patients with low antigen levels tend to have low post-treatment Tcon/Treg ratios (Figure 3D). Thus, we counted the post-treatment Tcon/Treg ratios in 23 SLE patients with different antigen titers. First, in the clinic, IL-2 therapy (1 million IU rhIL-2 every other day corresponding to 0.02 nM/h IL-2) reduces the Tcon/Treg ratios of 23 SLE patients, where the mean value of pre-treatment Tcon/Treg ratio is 8.43 and the mean value of post-treatment Tcon/Treg ratio is 4.78 (Figure 3E). The result suggests that 0.02 nM/h IL-2 is in the IDTW and can significantly reduce the Tcon/Treg ratio. Second, we found that the post-treatment Tcon/Treg ratios of SLE patients are correlated with the dsDNA titers (Figure 3F). The lower the dsDNA titer of an SLE patient, the lower the post-treatment Tcon/Treg ratio after IL-2 treatment, which is qualitatively consistent with our prediction (Figure 3D). Overall, our mathematical model is semi-quantitatively in agreement with clinical trials. Therefore, it is feasible to simulate the treatment process of SLE patients and give the most appropriate treatment regimen by our mathematical model. However, a more accurate mathematical model needs to be supported by larger-scale clinical data, and this is the direction of our future work.

Pre-treatment Treg of SLE patients can be a predictor of IL-2 treatment efficacy

Here, we discussed the heterogeneity of SLE patients. Clinical data suggest that the administration of IL-2 is not effective in all SLE patients (He et al., 2016). Failure of IL-2 therapy may be caused by administering IL-2 outside the IDTW, or results from individual factors, such as antigen level, gender, genetic factor, and age. Thus, it is important to judge whether the SLE patient is suitable for IL-2 treatment or not before treatment. Based on our model, we can predict an individual’s suitability for IL-2 therapy by classifying SLE patients before treatment.

To simulate heterogeneity among SLE patients, we randomly sampled 5,000 sets of all parameters in the range of one-fifth to five times the standard parameter set using Latin hypercube sampling. We simulated the post-treatment Tcon/Treg ratios to continuous IL-2 doses for each set of parameters. We found that these curves of Tcon/Treg ratios in response to IL-2 after PCA downscaling could be divided into two groups (k-means, Figure 4A). One group had a wide IDTW and a low post-treatment Tcon/Treg ratio. This was termed the IL-2-responsive group. Another group had a narrow IDTW and limited decrease in post-treatment Tcon/Treg ratio, or even no IDTW when the proliferation rate of Tcon was high and the proliferation rate of Treg was low, which is termed the IL-2-unresponsive group.

Figure 4

Classification of patients based on pre-treatment Treg

(A) Subjecting the response curve of Tcons of 5,000 parameter sets to PCA downscaling, these samples, each denoting a patient, can be divided into two categories (red and gray). Then, we labeled the PCA results by both the width of IDTW ( 0.05 nM/h, 0.05 nM/h, and no IDTW) and minimal post-treatment Tcon/Treg ratio ( 50% pre-treatment Tcon/Treg, pre-treatment Tcon/Treg and pre-treatment Tcon/Treg).

(B) Post-treatment Tcon/Treg ratio with exogenous IL-2 in the responsive group (red line) and unresponsive group (gray line) after Min-max normalization. The solid line represents the mean value of Tcon/Treg ratios, and the shaded part represents the standard deviation.

(C) The mean value and standard deviation of 10 parameters differed most between the responsive and unresponsive groups. These parameters are normalized by the mean values, respectively. The error bars denote the standard deviation.

(D) Mean proportion of pre-treatment Treg (Treg0) of the responsive and unresponsive groups (Mann–Whitney U test, p  0.001). The error bars denote the 95% confidence interval.

(E) The left panel shows the simulated result in which more samples in the Treg0 low group belong to the responsive group than those samples in the Treg0 high group. Clinically, the right panel shows that around 66.7% of patients in the Treg0 low group achieved SRI-4, but only 25.0% of patients in Treg0 high group achieved SRI-4 (Mann–Whitney U test, p = 0.044). The proportion of Treg0 less than 13.35% belongs to the Treg0 low group, and the proportion of Treg0 higher than 13.35% belongs to the Treg0 high group.

To investigate the main determinant in our patient classification system, we compared the parameters of these two groups and listed the ten most variable parameters between the two groups (Figure 4C). These parameters mainly influence pre-treatment Tcon and Treg density (named Tcon0 and Treg0) in SLE patients, including the differentiation rates of Tcon and Treg ( and ), proliferation rates of Tcon and Treg ( and ), the environmental carrying capacities of Tcon and Treg ( and ), the proliferation coefficients of Tcon and Treg induced by IL-2/IL-2R complex and average surface densities of IL-2R on Tcon and Treg . When combined with the fact that Treg0 in the IL-2-responsive group was lower than that in the IL-2-unresponsive group (Figure 4D), we inferred that Tcon0 and Treg0 may be good indicators to predict the classification of SLE patients before treatment.

We verified our speculations from both clinical data and simulations. In the clinical results, patients whose Treg0 in CD4+ T cells was less than 13.35% could be classified as the Treg0 low group, while patients whose Treg0 in CD4+ T cells was higher than 13.35% could be classified as the Treg0 high group. By analyzing the clinical data of 29 SLE patients (21 Treg0 low patients, 8 Treg0 high patients; see Tabel S1), we found that 66.7% of patients in the Treg0 low group could reach SLE Responder Index-4 (SRI-4), a measure of clinically meaningful improvement, after IL-2 treatment, while only 25.0% of patients in the Treg0 high group could reach SRI-4 (Figure 4E). In the simulation, we also showed that more patients with lower Treg0 belong to the IL-2-responsive group (Figure 4E). Thus, both theoretical prediction and experimental results suggest that patients with lower Treg0 are more likely to have a better therapeutic effect under IL-2 administration.

Collectively, these results indicate the feasibility of classifying SLE patients according to responsiveness to IL-2 exogenous treatment beforehand by measuring Treg0. In short, this means we can predict which SLE patients will respond well to IL-2 exogenous therapy and who will not. For those patients who are classified as unresponsive, traditional drug treatment, such as glucocorticoids, can be selected.

Two-variable model reveals the mechanism of treatment failure

The previous study indicated that high-dose IL-2 immunotherapy can lead to toxic adverse effects in the treatment of certain metastatic cancers (Krieg et al., 2010), herein we attempted to give a dynamic understanding of the mechanism of high-dose IL-2. As our SLE model contains many complex interactions that interfere us to perform more dynamic analysis of the IL-2 therapy, including nullcline analysis, etc. We developed a simplified model, containing only Tcon and Treg cells, using a modified Lotka–Volterra (LV) equation (Zhang et al., 2017), as follows:

To distinguish from the variables in the previous model, we used and to represent the density of Tcons and Tregs (The description can be found in STAR Methods). In the two-variable model, we focused on the competition between Tcon and Treg for IL-2 and the suppression of Tcon by Treg, which is the basis of the Tcon-Treg network. For simplicity, we did not consider the binding process of IL-2R to IL-2 and did not separate Treg into Treg and Treg according to their functions (Figure 5A).

Figure 5

Two-variable model reveals the mechanism underlying treatment failure with high-dose IL-2

(A) Schematic diagram of the two-variable model.

(B and C) shows the evolution of and (yellow line) under 0.01 (B) and 0.10 nM/h (C) IL-2. The gray dot denotes the pre-treatment state of the system, the yellow dot denotes the post-treatment state of the system under 0.01 nM/h IL-2, and the red dot denotes the post-treatment state of the system under 0.10 nM/h IL-2. The black line denotes the nullcline of , and the blue line denotes the nullcline of . For comparison, the nullclines of and without IL-2 are also shown by the black and blue dotted lines. The shadow part is the IDTW, which is defined by the 5.5 Tcon/Treg ratio.

(D) shows the post-treatment state of the system with increasing IL-2 doses. The gray dot denotes the pre-treatment state of the system, the yellow dot denotes the post-treatment state of the system under 0.01 nM/h IL-2, and the red dot denotes the post-treatment state of the system under 0.10 nM/h IL-2. The shadow part is the IDTW, which is defined by the 5.5 Tcon/Treg ratio.

(E) Dynamic landscape of treatment without IL-2, with low-dose IL-2 (0.01 nM/h), and with high-dose IL-2 (0.10 nM/h). The horizontal coordinates denote the density of and , and the vertical coordinate represents the doses of . The black lines with arrows denote the paths from the original attractor to low- and high-dose attractors after treatment. The color bar indicates the value of , where is the stationary probability distribution of the system and is used to define the landscape.

We demonstrated the evolution of Tcon and Treg cells on the phase space with 0.01 and 0.10 nM/h IL-2 (Figures 5B and 5C) and found that the directions of their evolutionary trajectories are different. The yellow point in Figure 5B denotes the steady-state with 0.01 nM/h IL-2, which is located in the IDTW, and the red point denotes the steady-state with 0.10 nM/h IL-2. The nullclines of and demonstrates the formation of the IDTW. When the IL-2 dose is low, the nullcline of is sensitive to IL-2, while the nullcline of is not. Thus, the nullcline of will move up, and the fixed point will evolve from the initial state to the yellow point in the IDTW. When administering high-dose IL-2, the nullcline of becomes sensitive to IL-2, while the nullcline of does not. Thus, the nullcline of will move right and push the fixed point outside the IDTW. Besides, we added Gaussian white noise to the variables and constructed the pseudo landscape in Figure 5E (Lv et al., 2015; Wang et al., 2008), which illustrates the stationary probability distribution of the system. High-dose (0.10 nM/h) IL-2 immediately causes the system to evolve from the initial point to the steady-state outside the IDTW without passing anywhere in the IDTW (Figures 5D and 5E), which means that high-dose administration shields the IDTW and does not help in the treatment of SLE.

Discussion

The immune system is a complex system involving different types of lymphocytes and cell-to-cell interactions with positive and negative feedback loops (Thurley et al., 2018). These factors determine the dynamic processes of the immune response and the pathogenesis of the immune disease. SLE, as a serious autoimmune disease, is characterized by the imbalance of regulatory and conventional T cells. The previous study has performed clinical trials to investigate the dose range of IL-2 therapy (Humrich et al., 2019). However, owing to ethical restrictions, the current IL-2 dose studies have the limitation of small samples, and high-dose IL-2 trials could not be performed in SLE patients either. Thus, the mathematical model can be a compliment to clinical trials. The quantitative study of SLE progression and its treatment will help in understanding the underlying mechanisms of the disease and improving treatment plans.

Herein, we built a semi-quantitative mathematical model to simulate the dynamic behavior of Tcon, Treg, and NK cells in SLE patients with IL-2 treatment. We identified the IDTW for SLE patients when treated with IL-2 and demonstrated the formation processes of IDTW. We found that patients with low antigen levels tend to have a wider range of IDTW and better therapeutic effects. Further, based on experimental results and our model, we divided patients into two groups, responsive and unresponsive, according to the response curves of Tcon with different IL-2 doses. Patients with low pre-treatment Tregs are more likely to be in the IL-2-responsive group, which is consistent with clinical findings.

We aimed to understand the mechanism of IL-2 therapy and improve the current treatment strategies. Although we proposed a model for the prediction and classification of patients based on pre-treatment Treg, this index ignores other immunological variations; therefore, its predictive accuracy still needs improvement. A possible solution could entail measuring specific parameters for each patient according to blood samples and modifying our model accordingly, which would provide for a more personalized approach. For example, we found that patients with high Tcon0 tend to have a higher likelihood of reaching SRI-4 (Figure S1). An indicator that includes both Tcon0 and Treg0 may be a better predictor of treatment effect. Ideally, combining clinical trials with the mathematical model might yield the best pathway to precision medicine.

Limitations of study

Our study does ignore the effects of the age, race, and gender on SLE disease. SLE mainly affects women in their 30s, which may be attributed to the effects of female sex hormones on the immune system. Age differences might influence the regulation naive T cell proliferation in the model. Race and gender differences, however, are difficult to build into our model, which will serve as a focus for subsequent work. Another drawback is our neglect of the IL-2R expression pathway. IL-2R on the surface of Treg cells is expressed continuously, but the IL-2R of Tcon will express only under the stimulation of IL-2. As we focused on the interaction between Tcon and Treg, we ignored the IL-2R expression pathway. Finally, a whole network of antigen presentation may consist of B, NK and other cells.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources should be directed to Fangting Li (e-mail: lft@pku.edu.cn).

Material availability

This study did not generate new unique reagents.

Method details

Clinical trial participants

This was a mathematical model study using published data from several previous clinical studies (NCT02465580, NCT02932137, NCT02084238 and DRKS00004858) of low-dose IL-2 in SLE patients. Detailed information on study designs, inclusion/exclusion criteria, patient demographic characteristics, studied drugs, concomitant medication and clinical/laboratory examinations for each completed study have previously been published (He et al., 2016, 2020; Humrich et al., 2019). These studies were conducted in accordance with the Declaration of Helsinki, the International Conference on Harmonization Guidelines for Good Clinical Practice, and were approved by ethics committees of the affiliations in which the studies were performed.

For the experiment of IL-2 production by immune cells (Figure 1C), the PBMCs from eight SLE patients and five healthy controls were obtained and tested by flow cytometry. Written informed consents were obtained from these participants. The experimental protocol followed the guidelines of the Declaration of Helsinki and was approved by the Human Ethics Committee of Peking University People’s Hospital (Beijing, China).

Flow cytometry analysis

Peripheral blood mononuclear cells (PBMCs) in SLE patients and healthy donors were analyzed by flow cytometry (FACSAria II; BD Biosciences). To analyze IL-2 production, PBMCs were stimulated and incubated for 4 h with 50 ng/mL phorbol myristate acetate, 1 g/mL ionomycin (both from Sigma-Aldrich), and 1 L/mL GolgiStop (BD Biosciences). After stimulation, the cell surface markers CD4, CD3, CD8, CD56, CD25 and CD127 were stained with the fluorescence labeled monoclonal antibodies. Next, the cells were fixed and permeabilised for 30 min at 4°C at dark with BD FACS fixation and permeabilisation buffer set (BD Biosciences). Then the cells were stained with PE labeled anti-human-IL-2 antibody. Proportion of IL-2 expressing cells in CD4+ effector T cells (defined as CD4+ CD8− CD127 CD25), double-negative T cells (DNT cell, CD3+ CD4− CD8−), NK cells (CD56+ CD3−) and NKT cells (CD56+ CD3+) was analyzed using FlowJo v10 software (Tree Star). Antibodies used for flow cytometry included APC-H7-anti-human CD3 (Biolegend), Percp-anti-human CD8 (Biolegend), FITC-anti-human CD4 (Biolegend), BV421-anti-human CD25 (Biolegend), BV605-anti-human CD127 (Biolegend), APC-anti-human CD56 (Biolegend) and PE-anti-human IL-2 (Biolegend).

Volume effect factor of IL-2 releasing

Here we derived the volume effect factor as follows. When the transmembrane transport of IL-2 is not considered, in a very short time d, the increment of extracellular IL-2 molecules in is equal to the increment of IL-2 molecules produced by all Tcon cells , . We divided both sides by , and the left side is the incremental concentration of extracellular IL-2 (d ), which is given by d . Since represents the incremental concentration of IL-2 produced by all Tcons, and it can be written as . Note that  = ( denotes the volume of one single cell) and that the term of IL-2 generated by all Tcons can be written as:

We simplified to in our model, the term describes IL-2 secreted by Tcons can be written as . Because of the low cell density in peripheral blood, is approximated as 1, thus approximated to mm3. In the vicinity of tissues or organs, however, the value of is greater than 1. We discussed the effects of environments on SLE pathogenesis in Figure S9. Similarly, the volume effect factor exists on every term involving the release from intracellular to extracellular.

Conversion factor of IL-2R

For simplicity, we assumed that IL-2R and IL-2/IL-2R complex are evenly distributed intracellular through introducing the conversion factor . Since this assumption involves changes in the unit of IL-2R, we introduced the notation to indicate the value of the variable at a given unit. Within a single cell volume (m3), we set the surface density of IL-2R as (m−2) and the concentration of IL-2R in a single cell after conversion is (nM). The number of IL-2R in a single cell () can be derived by the concentration of IL-2R as , where is the Avogadro constant. Also, the number of IL-2R can be derived by the surface density of IL-2R as , where (m) is the radius of a single cell. Based on a conservation of the number of IL-2R on a single cell, we had the equation . Here, we defined as the conversion factor ( is approximately equal to 1 m2 nM when setting to 5 m), and the equation above can be simplified as

Molecular perspective of IL-2 binding process

We deduced the process of IL-2 binding to IL-2R at the molecular level. Total cells are distributed evenly in the system, where indicates the total volume in peripheral blood, and indicates the volume of one single cell with radius . We considered the probability of a single IL-2 molecule binding with the IL-2R on a single cell, and here we viewed IL-2 approximately at a particle point. The volume of free space in peripheral blood can be written as

In the stationary reference frame relative to IL-2, each single cell moves at a relative speed . In interval , the volume covered by the movement of cell is

IL-2 molecules can appear at any place in free space with the same probability. When IL-2 overlaps with space covered by T cell in , we called the IL-2 molecule collided with the cell. The probability that an IL-2 molecule collides with one cell in d time is written as:

It should be emphasized that we default to the cell not colliding with other cells in time. We can always find the interval for which adjacent cells will not collide with each other within a limited time, as long as 0. Assuming that the relative speed of IL-2 molecule to cell is , the number of collisions of all IL-2 molecules with one single cell in time iswhere is the number of IL-2 molecules, is the average relative speed of all IL-2 to cell . It can be considered that each cell has the same value of the average relative speed of all IL-2 molecules, thus we used to replace . Notice that , we can simplify the equation above and get the average times of collision as:

We have derived the number of collisions between one single cell and IL-2 molecules per unit of time. Let us denote the reaction cross-section of an IL-2 molecule and an IL-2R to , and the number of IL-2R molecules per unit surface area on a cell to , where is the number of IL-2R on a single cell. Assuming IL-2 molecules will appear anywhere on the surface with the same probability, and the probability of a single IL-2 molecule binding to IL-2R among one collision is

Finally, the times that IL-2 binding with IL-2R on one Tcon cell per unit of time can be written as:where denotes the binding rate of IL-2 and equals to .

Rapid equilibrium assumption of sCD25/IL-2 complex

sCD25 is the chain of IL-2R, which is widely distributed in the blood of SLE patients. We set the binding rate of sCD25 and IL-2 as , the dissociation rate of sCD25/IL-2 complex as and the degradation rate of sCD25/IL-2 complex as . The reaction of the binding process of sCD25 and IL-2 can be represented by

Assuming that the binding process of IL-2 and CD25 is faster than the degradation process. The sCD25/IL-2 concentration can be described by the ordinary differential equation as

In the steady-state, we have the equation

We assumed that sCD25 is conserved in the system, i.e., . Thus, we can derive the binding rate of sCD25 and IL-2 aswhere . Setting to 0.83 h−1, to 1.8 h−1 and to 37.2 h−1 (Busse et al., 2010), approximates to 0.07 nM. Similarly, the dissociation rate of sCD25/IL-2 complex is . Considering the dissociation of sCD25/IL-2 complex, the decrement of IL-2 after correction can be written as

The decrement of IL-2 after correction is around 0.672 times the binding rate of sCD25 and IL-2. Therefore, considering the dissociation of the sCD25/IL-2 complex would bring a 0.672-fold correction to the average concentration of sCD25, which is easy to simulate by adjusting the parameter value. Besides, the mechanism of sCD25 production is not fully understood and the interactions between sCD25 and immune cell subsets are complex, which needs further study. Also, to avoid introducing too many variables, we ignored the influence of IL-2 on the density of sCD25 and take the concentration of sCD25 as a parameter in our model rather than a variable to represent the competitive effect of extracellular sCD25 on IL-2.

Parameter setting and estimation

Then we demonstrated the estimation of these parameters.

Part 1 estimates the parameters in cell proliferation. Under antigen stimulation, the apoptosis rate of CD4+ T cells vary from to day−1 and the apoptosis rate of APCs vary from to day−1 (Price et al., 2015), we chose 0.01 h−1 as the apoptosis rate of CD4+ T cells, APCs and NK cells ( and ). Based on the facts that the average number of CD4+ T cells in blood is around 1000 mm−3 (Shete et al., 2010) and the range of TNs is 216-1059 L−1 blood (Bains et al., 2009), we calculated the renewal rate of TNs (), the differentiation rate ( and ) through three assumptions. First, in SLE patients, Treg is dysfunctional and the suppression of Treg to Tcon is negligible. Second, in the case of IL-2 deficiency, the proliferation of Tcon induced by the IL-2/IL-2R complex is negligible. Third, under adequate antigen stimulation, the Hill functions induced by antigens are close to 1.

We derived the differentiation rate of TNs as follows. Setting d , the steady-state of Tcons can be written as according to our assumptions. Setting to 1000 mm−3 (Shete et al., 2010) and to 0.01 h−1, it satisfies that . Since the range of TN is 216-1059 L−1 blood (Bains et al., 2009), the differentiation rate of TNs to Tcon () is approximately 0.01–0.05 h−1 and we set the value as 0.018 h−1 in our model. Considering that Treg density is roughly one-10th of CD4+ T cells in SLE patients before treatment (He et al., 2016), we estimated the differentiation rate of TNs to Treg () as /9 = 0.002 h−1.

Then, we derived the renewal rate of TNs. Setting d , under adequate antigen stimulation, the steady-state of TNs satisfy , i.e., . Thus, the range of renewal rate of TN is 6.48–31.77 mm3 h−1 and we set the value as 20 mm3 h−1 in our model. Since the pretreatment NK cell density is around 60 mm−3 (Humrich et al., 2019), we calculated the renewal rate of NK cell () as mm3 h−1. Considering the T/APC ratio is around 10:1 (Friedl and Gunzer, 2001), the steady-state of APCs is around 100 mm−3 and we estimated the renewal rate of APC () to  = 1 mm3 h−1. The range of T cell proliferation is 0-6 days−1 (Burroughs et al., 2006), thus, we estimated the proliferation rate of APC cell (), Tcon (), Treg (, ) and NK cell () both to 0.05 h−1. We estimated the environmental carrying capacities of Tcons, Tregs, NK cells and APC cells (, , and ) based on the previous work (Shete et al., 2010; He et al., 2016; Humrich et al., 2019; Friedl and Gunzer, 2001), and set to 3000 mm−3, to 400 mm−3, to 400 mm−3 and to 200 mm−3. Other parameters related to cell proliferation (, , and ) were estimated.

Part 2 estimates the parameters related to IL-2 secretion. In a healthy person, the secretion rate of IL-2 is 0-22000 molecules∗cell−1∗h−1 (Busse et al., 2010) and we set the value to 10,000 molecules∗cell−1∗h−1. In a system with a volume of 1 L, the total number of cells in our model is . The increase in IL-2 concentration per unit time in the system can be written as nM, where is the Avogadro constant. Thus intracellularly, the expression rate of IL-2 of Tcons in healthy person satisfies according to the conservation of IL-2 molecular number. Thus, can be simplified as , which approximately equals 33 nM/h (0–72.6 nM/h). Considering the secretion capacity of SLE patients are weaker than healthy people (Figure 1C), so we estimated the IL-2 secretion rate of SLE patients ( and ) to 1 nM/h. The degradation rate of IL-2 is obtained from Busse et al.’s work (Busse et al., 2010).

Part 3 estimates the parameters related to IL-2R as follows. The binding rate , dissociation rate of IL-2/IL-2R , and degradation rate of IL-2/IL-2R complex were obtained from Busse et al.’s work (Busse et al., 2010). Considering the IL-2R on NK cells is IL-2R , we assumed the binding rate of IL-2R on NK cells is one-third of to differ with IL-2R on Tcons and Tregs. Besides, since sCD25 is the alpha chain of IL-2R, we estimated the binding rate of IL-2 to sCD25 is one-third of . Previous work gave the number of IL-2R on Tregs as 104 /cell (Busse et al., 2010). Considering that Tregs are more competitive for IL-2 than Tcons (Abbas et al., 2018) and the density of Tregs is only one-10th of CD4+ T cells (He et al., 2016), we estimated that the IL-2R on Tregs are two orders higher than IL-2R on Tcons. Thus, we set the average number of IL-2R on one Treg cell to 104/cell, the average number of IL-2R on one Tcon and NK cell to 102/cell (We investigated the influence of the choice of the average number of IL-2R on a cell in Figure S5). The cell surface area is about 300 m2 (Busse et al., 2010), thus the average density of IL-2R on Tcons () and NK cells () are 0.3 m−2, and the average density of IL-2R on Tregs () is 30 m−2 sCD25 is also present in the serum of SLE patients, at a concentration of 271.4 pg/mL for SLE patients (El Shafey et al., 2008). Due to sCD25 is a 24.8 kDa protein (http://tds.tonbobio.com/tds-21-9241.pdf), we estimated the concentration of sCD25 () at 0.01 nM. Last, the proliferation coefficients induced by IL-2/IL-2R are estimated based on the average density of IL-2R.

The level of dsDNA can roughly represent the level of self-antigen. We have shown the titer of dsDNA in SLE patients in Figure 3. It is reasonable to think self-antigen in SLE patients is saturated and enough to activate the immune system. To facilitate the calculation, we set the patient’s self-antigen to 5 A.U. and the antigen-induced proliferation coefficient to 1 A.U.

We used recombinant human IL-2 (rhIL-2) in the clinic to treat SLE at doses ranging from 0.2 million IU to 2 million IU every two days in the past research (He et al., 2016). The biological activity of rhIL-2 is around 107 IU/mg, and the molecular weight of rhIL-2 is 15 kDa (https://www.rndsystems.com/cn/products/recombinant-human-il-2-protein_202-il). Assuming that the total volume of human blood is 5 L, the concentration of 1 million IU of rhIL-2 in the blood after complete absorption is 1 nM. We usually used rhIL-2 every two days, so that the flux of exogenous IL-2 equals 0.02 nMh−1. The clinical dose corresponds to a flux ranging from 0.004 to 0.04 nMh−1.

Two-variable model

Interactions among Tregs, Tcons and exogenous IL-2 are the central processes during SLE treatment. To further investigate the dynamic behavior of IL-2 treatment, we constructed a simplified two-variable model including Tcons and Tregs. We developed the following mathematical equations using a modified Lotka-Volterra (LV) equation (Zhang et al., 2017), as follows:

To distinguish from the variables in the previous model, we used and to represent the density of Tcons and Tregs here. The first term in Equation 11 describes the renewal of Tcons, and the second term describes the apoptosis of Tcons. The third term is the LV equation to show the competition between Tcons and Tregs for IL-2, and the last term shows the suppression of Tcons by Tregs. Equation 12 gives the density of Tregs, where the first term shows the renewal of Treg and the second term shows the apoptosis of Tregs. Similar to Tcons, we used the LV equation to describe the competition for IL-2 between Tcons and Tregs.

Part of the parameters have the same meaning as in Table 1, thus we used the same values as in Table 1, i.e., the degradation rate of Tcons  = 0.01 h−1, the degradation rate of Tregs  = 0.01 h−1, the environmental carrying capacity of Tcons  = 3000 mm−3, the environmental carrying capacity of Tregs  = 400 mm−3, the suppression rate of Tcon  = 0.02 h−1 and the suppression coefficient of Tcons induced by Treg  = 500 mm−3. Besides, we set the renewal rate of Tregs  = 1 mm−3 h−1 to make the pretreatment Treg value  = 100 mm−3 h−1 (He et al., 2016). Also, we set the renewal rate of Tcons  = 10 mm−3 h−1 to make the pretreatment Treg ratio to 10.4% (He et al., 2016). We assumed that the proliferation rate of Tcon depends on the dose of exogenous IL-2 in the form of , where  = 0.05 h−1 (the same as Table 1) and = 0.05 nMh−1 denotes the proliferation coefficient of Tcons induced by IL-2. Similarly, the proliferation rate of Treg depends on the dose of exogenous IL-2 in the form of , where  = 0.05 h−1 (the same as Table 1) and = 0.01 nMh−1 denotes the proliferation coefficient of Tregs induced by IL-2. The matrix here indicates the competitive interactions between Tcons and Tregs, whereas and indicate the effect of Tcons and Tregs on themselves, and we set and equal to 1. Furthermore, the competition interactions of Treg to Tcons were relatively strong, and we estimated  = 5 and  = 0.1.

Quantification and statistical analysis

Data in Figure 1C were expressed as the median and range for non-normally distributed data. The Mann-Whitney U test was performed for non-parametric data.

Data in Figure 3F were expressed as median with IQR. The paired t-test was performed.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
APC-H7-anti-human CD3	BD Biosciences	Cat. No. 560176; RRID: AB_1645475.
Percp-anti-human CD8	Biolegend	Cat. No. 344707; RRID: AB_1967122.
FITC-anti-human CD4	Biolegend	Cat. No. 317408; RRID: AB_571951.
BV421-anti-human CD25	Biolegend	Cat. No. 302630; RRID: AB_11126749.
BV605-anti-human CD127	Biolegend	Cat. No. 351334; RRID: AB_2562022.
APC-anti-human CD56	Biolegend	Cat. No. 362504; RRID: AB_2563913.
PE-anti-human IL-2	Biolegend	Cat. No. 500307; RRID: AB_315094.
Chemicals, peptides, and recombinant proteins
phorbol myristate acetate	Sigma- Aldrich	P1585
Ionomycin	Sigma- Aldrich	I3909
GolgiStop	BD Biosciences	554724
Critical commercial assays
BD FACS fixation and permeabilization buffer set	BD Biosciences	Cat. No. 554714; RRID: AB_2869008.
Deposited data
Tcon and Treg data	(He et al., 2016)	https://doi.org/10.1038/nm.4148
NK cell ratio	(He et al., 2020)	https://doi.org/10.1136/annrheumdis-2019-215396
NK cell number	(Humrich et al., 2019)	https://doi.org/10.1016/S2665-9913(19)30018-9
Software and algorithms
SLE model	This paper	https://doi.org/10.5281/zenodo.6922465