Predictive Models of Type 1 Diabetes Progression: Understanding T-Cell Cycles and Their Implications on Autoantibody Release

Defining the role of T-cell avidity and killing efficacy in forming immunological response(s), leading to relapse-remission and autoantibody release in autoimmune type 1 diabetes (T1D), remains incompletely understood. Using competition-based population models of T- and B-cells, we provide a predictive tool to determine how these two parametric quantities, namely, avidity and killing efficacy, affect disease outcomes. We show that, in the presence of T-cell competition, successive waves along with cyclic fluctuations in the number of T-cells are exhibited by the model, with the former induced by transient bistability and the latter by transient periodic orbits. We hypothesize that these two immunological processes are responsible for making T1D a relapsing-remitting disease within prolonged but limited durations. The period and the number of peaks of these two processes differ, making them potential candidates to determine how plausible waves and cyclic fluctuations are in producing such effects. By assuming that T-cell and B-cell avidities are correlated, we demonstrate that autoantibodies associated with the higher avidity T-cell clones are first to be detected, and they reach their detectability level faster than those associated with the low avidity clones, independent of what T-cell killing efficacies are. Such outcomes are consistent with experimental observations in humans and they provide a rationale for observing rapid and slow progressors of T1D in high risk subjects. Our analysis of the models also reveals that it is possible to improve disease outcomes by unexpectedly increasing the avidity of certain subclones of T-cells. The decline in the number of -cells in these cases still occurs, but it terminates early, leaving sufficient number of functioning -cells in operation and the affected individual asymptomatic. These results indicate that the models presented here are of clinical relevance because of their potential use in developing predictive algorithms of rapid and slow progression to clinical T1D.

S1.1. ONE-CLONE MODEL S1.1.1. Model Rescaling. We first start by nondimensionalizing the model (4)- (9) in the main text in order to reduce the number of parameters and simplify the analysis.
Since T-cell activation is regulated by the level of pMHC expression level on APCs, we expect that dF dP ≥ 0.
Based on T-cell dose-response curves to growing level of pMHCs on APCs shown in [1], one can typically describe the function F(P) by a Hill function with a Hill coefficient n = 1, given by F (P) = P P + k , (S1.1) where k represents the expression level of pMHCs for 50% maximum activation of T-cells. The larger the value of k, the higher the expression level of pMHCs required for T-cell activation, which implies that k −1 can be used to quantify T-cell avidity. The empirical determination of the Hill function that approximates T-cell activation (by measuring INF γ secretion from CD4 + T-cells in response to variations to GAD-expression level on APCs) was also done in [2]. As explained in [3], the value of k for GAD-reactive T-cells was found to lie between [0.02,0.18] µM for highavidity T-cells [2], but this range was then extended to [0. 1,9] µM for low-avidity T cells. We take advantage of a broad spectrum for k (∈ [0, 20] µM) in our analysis of the model (4)- (9) to study the effects of T-cell avidity on disease onset and progression.
The function H 1 in (8) determines β-cell killing by (CD8 + and CD4 + ) T-cells occurring at an effective rate κ in the range [10 −11 , 10 −7 ] (day·cell) −1 (which is expected to remain approximately constant throughout disease progression for a given individual [4,5]). On the other hand, the function H 2 in (9), which describes pMHC production from processing β-cell specific proteins in APCs, is assumed to be proportional to H 1 (T, β). In other words, H 2 (T, β) = R H 1 (T, β), (S1. 2) where R is the production rate of pMHC per T-cell per β-cell. By using mass-action kinetics to describe the autoimmune attack by CD8 + and CD4 + T-cells and by assuming that β-cell killing is 1 saturated by β-cell number, we arrive at the following formalism for H 1 where µ is the saturation parameter per β-cell for β-cell killing. In some special cases, we may assume that the saturation parameter is identically zero to make H 1 a bilinear function of T-cell and β-cell population sizes. β-cell renewal (via β-cell replication or neogenesis) can be described by (see [4]) where k β denotes the number of β-cells required for 50% maximal renewal. In Eqn. (4), the source term Σ T represents the input of naïve cells from the thymus, which does not necessarily remain constant throughout disease progression. The assumption here is to take the source term Σ T to satisfy Σ T = σ F (P). (S1.5) By making use of Eqns. (S1.1)-(S1.5), we can give a detailed and explicit mathematical description to the one-clone model (4)-(9) as follows dT where s is the maximal rate of β-cell renewal per day, whereas δ β is the turnover (mortality) rate of β-cells per day. The new model (S1.6)-(S1.11) differ from the one presented in [3] in three main aspects; namely, the presence of a term describing thymus input Σ T , the dependence of pMHC production rate on κ (T-cell killing efficacy), and the dependence of B-cell activation on the pMHC-dependent Hill function F . These modifications make the current model more physiological, especially the latter one which makes T-and B-cell avidities more closely correlated via the parameter k than previously assumed.
In order to further reduce the complexity of the model or extent of the analysis, we simplify the model (S1.6)-(S1.11) by making the following substitutions where β 0 is the initial number of β-cells right before the autoimmune attack (with the ratio β/β 0 denoted by a non-italicized variable β) and For simplicity, we assume that s = µ = δ β = 0. By applying the substitutions introduced in (S1.12) to the model (S1.6)-(S1.11) and using non-italicized font for the newly generated parametric quantities hereafter, we obtain where the new parametric quantities are given by Thus the number of parameters that need to be estimated have been reduced by introducing these parametric quantities. S1.1.2. Parameter Estimation. For the (scaled) one-clone model (S1.13)-(S1.18), we adopt the parameter values that appeared in [3][4][5][6][7][8][9][10] to perform our simulations. These parameters were estimated using a combination of fit and parsimony. For the definitions, values and units of all these parameters, please see Table S1.1.
In brief, the kinetic parameters associated with B-cell and plasma-cell expansion and turnover were previously determined using in vitro data from [10], whereas those associated with autoantibody release and degradation were determined based on IgG antibody isotype data from [11] and model rescaling (see [3] for more details).  S1.1. Parameter values of the (scaled) one-clone model (S1.13)-(S1.18).
As for the kinetic parameters of T-cells, β-cells and autoantigen(s), we rely on previously estimated values from [3][4][5][6][7][8]. These estimates were generated using experimental data available in the literature on T1D. Where relevant, we assign similar values to the parameters listed here.
Because of heterogeneity in the avidities and killing efficacies of the different clones of T1Dspecific T-cells (and B-cells), physiologically reasonable ranges for these parameters are instead considered in the simulations shown here. Experimental data from [1,2] were used to determine such ranges. S1.1.3. Phase-Plane Analysis. In this section, we will study the global dynamics of the one clone model (4)-(9), described in the paper, using steady state and stability analysis. We will show that the full one clone model possesses a stable disease-free (healthy) steady state S 0 = (0, 1, 0, a, β ∞ , 0), with 0.3 ≤ β ∞ ≤ 1, that corresponds to healthy individuals, and a transient (quasi-stable) autoimmune state S 2 (with elevated levels of autoreactive T-cells, B-cells, plasmacells and autoantibodies) that corresponds to type 1 diabetic patients, both separated by a transient interior saddle point S 1 . Such behavior is similar to that observed in the one clone model of [3]. In the following, we use a reductionist approach to determine these steady states and analyze their stability properties.
We begin first by noting that due to homeostatic mechanisms, β-cell loss happens at a very slow time scale relative to T-cell dynamics and pMHC processing. Therefore, we may assume that the scaled variable β is roughly a constant, i.e., β = 1. Furthermore, since T and P are decoupled from B, P c and I g , we can ignore the differential equations of these latter variables and only focus on the reduced (scaled one-clone) two-variable subsystem To investigate the steady state dynamics of the reduced subsystem (S1.19)-(S1.20), we utilize the phase portrait in Fig. S1.1 showing the T -and P-nullclines defined by the equations T ′ = 0 and P ′ = 0, respectively. These nullclines subdivide the positive quadrant of the P, T -plane into regions where P ′ and T ′ have various signs. Figure S1.1 also shows the vector-field to illustrate the overall flow pattern of solution trajectories of the reduced model. The steady states of the reduced model lies at the intersection of these nullclines. The origin E * 0 = (0, 0) is a steady state representing the disease-free state, whereas the other steady states E * = (T * , P * ), if they exist, must satisfy P * = T * (S1.21) and where T * 0. Figure S1.1 shows that when this equation is satisfied, at most two more interior equilibria E * 1 = (T * 1 , P * 1 ) and E * 2 = (T * 2 , P * 2 ) must also be present, with the latter representing the autoimmune state. In fact decreasing the value of k gradually shifts the T -nullcline upward and increases the number of intersections between the two nullclines from zero to two. This is illustrated in Fig. S1.1 in which we see either no intersection for k > k c (gray dashed-dotted line), one intersection at k = k c (gray dashed line) and two intersections for k < k c (black solid line), where k c is the critical value of k at which the two points E * 1 and the autoimmune state E * 2 merge. The value of k c can be determined by solving for T * from (S1.22) as follows where ℓ = √ α − √ δ T (taken to be positive). Therefore to obtain one root to Eqn. (S1.23), k c must satisfy the quadratic equation (S1.24) If we ignore thymus input by setting σ = 0, we obtain k c = 1 (or ] Because the thymus input does not alter the dynamics of the model significantly [7], we only focus our analysis of the stability properties of all possible steady states hereafter and throughout the paper to the case when σ = 0. This is done by evaluating the Jacobian matrix of the reduced The phase portrait of the reduced (scaled) one-clone model (S1.19)-(S1.20), with σ = 0 and ℓ > 0, displaying the nullclines at different values of k.
The T -and P-nullclines are shown as solid black and gray lines, respectively, for k < k c = 1 (the P-axis is also a T -nullcline). Increasing the value of k shifts the Tnullcline downward, becoming first tangential to the P-nullcline at k = k c (dashed gray line), then lying entirely below the P-nullcline when k > k c (dashed-dotted gray line). The steady states, lying at the intersection of the T -and P-nullclines, E * 0 and E * 2 (stable nodes) are shown as black dots and E * 1 (saddle) as hollow squares. The stable manifold of the saddle point E * 1 (the black and gray dotted curves) is the separatrix that demarcates the basin of attraction of E * 0 from that of E * 2 .
(scaled) one-clone model (S1.19)-(S1.20) at these steady states, as follows The eigenvalues of the Jacobian matrix at the disease-free (healthy) state E * 0 = (0, 0) are both negative (λ 1 = −δ T and λ 2 = −δ P ), which means that E * 0 is always a stable node. The other two steady states E * 1 and the autoimmune state E * 2 , on the other hand, will only exist if k < k c = 1 and ℓ > 0. Assuming that these two conditions are satisfied and noting that the trace of the Jacobian matrix is always negative, it suffices to show that det(J(E * 1 )) < 0 and det(J(E * 2 )) > 0 to prove that E * 1 is a saddle and E * 2 is a stable node [12]. To achieve this goal, we evaluate the determinant of the Jacobian matrix and determine when it changes sign. Since every non-zero steady state (i.e., those that are not E * 0 = (0, 0), when σ = 0) of the reduced (scaled) one-clone model (S1.19)-(S1.20) must satisfy (S1.21) and (S1.22), then according to (S1.25), we have Thus det(J(T , P)) = 0, . This means that if T > T , as is the case for the autoimmune state E * 2 , then det(J(T , P)) > 0 and E * 2 is a stable node, whereas if T < T , as is the case for E * 1 , then det(J(T , P)) < 0 and E * 1 is a saddle. At T = T = T * (k c = 1), we have a saddle-node bifurcation as demonstrated in Fig. S1.1.
These latter stability results illustrate the local behavior of the reduced (scaled) one-clone model (S1.19)-(S1.20) near its steady states, when σ = 0. It remains to check if these results are globally valid. In the following analysis, we show that they are indeed global. This is done by establishing first that the solutions of the reduced (scaled) one-clone model, whose initial conditions belong to the space R + × R ≥0 , are bounded and nonnegative, and by taking advantage of Dulac's Criterion to rule out the existence of closed orbits in that space.
We now show that no periodic orbits exist using Dulac's Criterion.
Theorem S.2 [Dulac's Criterion, see [13], pp 202] Letẋ = f(x) be a continuously differentiable vector field defined on a simply connected subset R of the plane. If there exists a continuously differentiable, real-valued function g(x) such that ∇ · g(x) has one sign throughout R, then there are no closed orbits lying entirely in R.
By choosing the continuously differentiable vector field g(T , P) = 1/T in R + × R ≥0 , and applying Dulac's Criterion onẋ = (Ṫ ,Ṗ ) that satisfies the reduced (scaled) one-clone model (S1.19)-(S1.20) (with σ = 0 and ℓ > 0), we obtain The control condition ∇ · (gẋ) is always negative for T > 0 and P ≥ 0. As a result, there are no closed orbits lying entirely in the space (T , P) ∈ R + × R ≥0 . Furthermore, we already know from the previous result that solutions starting from R + × R ≥0 will never leave this space. Thus no periodic orbits would lie partially within R + × R ≥0 . These conclusions guarantee that the local stability properties of the steady states E * i , i = 0, 1, 2, are also global. [Similar result can be obtained for σ 0.] In this case, when k < k c , we have a bistable system with the stable manifold of the saddle point E * 1 (plotted as a black dotted line in Fig. S1.1) is a separatrix that divides the two basins of attraction of E * 0 and E * 2 . This configuration remains almost the same when k = k c . Here the separatrix (plotted as a gray dotted line) becomes the boundary between the basin of attraction of E * 0 and the basin of attraction of E * 1 = E * 2 , which is a half-stable steady state. It is important to point out that when considering the full model (S1.13)-(S1.18), the steady states E * i , i = 0, 1, 2, automatically translate into the steady states S i introduced at the beginning of this section. However, due to the fact thatβ < 0 for T > 0 in Eqn. (S1.17), we expect β(t) to be an exponentially decaying function of time, and that the two steady states S 1 and S 2 , which both possess elevated levels of T-cells, to be transient steady states, whenever k < k c . In other words, solution trajectories that start from the basin of attraction of S 2 will initially approach S 2 , but eventually turn towards S 0 when the slope of the hyperplane T = P/β becomes large enough that the two steady states S 1 and S 2 merge and disappear (i.e., k becomes larger than k c ). Similar behavior was previously observed in [3]. To illustrate the dynamics of the full (scaled) one-clone model in the presence of all of its components, we simulate system (S1.13)-(S1.18) in response to variations to two key parameters in the model; namely, the reciprocal of T-cell avidity k and T-cell killing efficacy κ (other parameter values are available in Table S1.1). In these simulations, the scaled variables are used as representatives to the original variables of system (S1.6)-(S1.11). Furthermore, because thymus input is ignored, the initial T-cell level T is taken to be nonzero but ≪ 1. Our goal in these simulations is to show that modifying the (scaled) one-clone model, by making B-cell activation depend on k (or k in the nonscaled model (S1.6)-(S1.11)) via the Hill function F(P) = P n P n + k n , (S1.27) where n = 1 is the Hill coefficient, will produce similar results to those obtained from the oneclone model analyzed in [3] and to uncover other aspects of the model. Figure S1.2(A1) shows the heat-map of the (scaled) level of β-cells (β) at steady state (i.e., after 30 years of the autoimmune attack) when k and κ are varied within the ranges [10 −4 , 1] and [10 −11 , 10 −7 ] (day·cell) −1 , respectively. We see that the magnitude of β-cell loss (shown as a gradual change in color from red to blue) increases steadily by increasing T-cell avidity and/or its killing efficacy in a manner identical to what have been observed previously in [3]. The presence of the red band on top of each panel demonstrates that if T-cell avidity is too small, then β-cells are safe from T-cell destruction regardless of the level of T-cell killing efficacy κ. Increasing the Hill-coefficient to n = 2, in panel (A2), or n = 3, in panel (A3), does not alter the heat-map of β significantly. The only noticeable difference we observe is the increase in the width of the red band as n increases. The increase in the steepness of the Hill function, described in (S1.27), for larger n means higher pMHC expression level on APCs is required for T-cell activation and thus wider red bands. The sudden change from red to blue when moving vertically downward across these red bands, however, is caused by the presence of bistable nodes whose basins of attractions are separated by the the stable manifold of the saddle point (the separatrix discussed in the previous section) which generates the boundary of these red bands.
We also observe on the left side of panels (A1-A3) red bands that are less darker than those on top. These red bands demonstrate that most β-cells survive the autoimmune attack whenever T-cell killing efficacy is very small (i.e., when κ ∈ [10 −11 , 5 × 10 −10 ] (day·cell) −1 ). For larger values of κ, on the other hand, a more significant decline in β is detected. In fact, increasing the value of κ beyond the critical threshold, highlighted by the thick black lines, labeled 0.3, makes β-cell survival below 30%, indicating clinical manifestation of the disease (through the appearance of T1D-associated symptoms). In the extreme cases identified in the bottom right corner of panels (A1-A3), the survival rate of β-cells is very small due to the effective autoimmune assault dominated by high-avidity, high-killing efficacy T-cells.
Given that B-cell activation also depends on the Hill function described by (S1.27) with n = 1, we further characterize the effect of k and κ on the level and survival of circulating autoantibodies. To do so, we plot in panels (B1-B3) the heat-map of I g at three successive time points: six months after the inception of the autoimmune attack (B1), at the clinical onset of the disease when the 0.3-threshold is crossed (which applies only to the points to the right of the black line) (B2) and at steady state (30 years after the inception of the autoimmune attack) (B3). Notice that model outcomes here are almost identical to those observed by the one clone model described in [3]. In brief, we observe four possible scenarios: I g becomes elevated without reaching diagnostic T1D (to the left of the black line), I g becomes elevated while reaching diagnostic T1D (to the right of the black line and with κ ≤ 10 −9 ), I g remains elevated until disease onset (to the right of the black line with κ ∈ [10 −9 , 8 × 10 −8 ]), I g remains undetectable throughout disease progression (the blue regimes close to the right edge of each panel). The high κ value in the latter case causes β-cell destruction to be too fast for autoantibody accumulation. Most of these outcomes are consistent with experimental observations in humans and animal models that screened positive for T1Dspecific autoantibodies [14,15].
To determine how fast the level of autoantibodies rise with respect to time while varying k and κ, we quantify in panels (B1-B3) the duration (in days) for the autoantibodies to reach the following detectability levels: 0.15 (C1), 0.55 (C2) and 0.95 (C3). As demonstrated by these panels and the color bar on the right, the rise in the level of autoantibodies in most cases is very fast and reaches its maximal level of 0.95 in less than 200 days after the engagement of T-cells in the destruction of β-cells. It should be mentioned here that the blue regime in panels (C2) and (C3) for large κ, indicates that the maximal detectability levels chosen (0.55 and 0.95, respectively) are never attained in these two cases. This result could be used as a criterion to determine when high risk subjects could be tested for autoantibodies as a diagnostic tool.
The examination of the levels of β-cells (β) and autoantibodies (I g ) in Fig. S1.2 can be better understood by tracking the time evolution of these two variables as well as T-cells (T ) in response to variations to the same key parameters k and κ. As shown in Fig. S1.3, the 30-year time evolution of T (A1-A4), I g (B1-B4) and β (C1-C4), after the start of the autoimmune attack, are plotted as heat-maps. These heat-maps are generated by taking k ∈ [10 −4 , 1] and choosing the following values for the killing efficacy κ: 10 −11 (A1-C1), 10 −10 (A2-C2), 10 −9 (A3-C3), 10 −8 (A4-C4) (day·cell) −1 . The one consistent feature observed across all these panels is the absence of any autoimmune attack at all time when k is close to 1 (i.e., T-cell avidity is too small to invoke an autoimmune response or even illicit T-cell expansion). However, when k is large enough, T-cell expansion becomes prominent as demonstrated by the fast increase (within few months) in T and the appearance of red regimes in panels (A1-A4). The increase in T does not always imply a significant loss in β-cells. In fact, panel (C1) shows that even after 30-years of follow up, we do not see much loss in β-cells, because κ is too small for T-cells to cause any harm. By increasing κ, β-cell loss becomes significant as demonstrated by the appearance of blue regimes in panels (C2-C4). The bigger the value of κ, the faster the loss in β-cells and the quicker the manifestation of the disease. T-cell survival, in these cases, is not maintained due to the decline in β-cell specific peptides that are produced from apoptotic β-cells required for T-cell activation. The rise and decline in T create these (red) "waves" that are induced by the transient bistability of the model. This explains the appearance of blue regimes in panels (A3-A4) in later years which are compatible with the blue regimes in panels (C3-C4). As for the level of circulating autoantibodies, panels (B1-B4) show that the time evolution of I g is similar to that of T , except for the delay in the rise of I g to its peak when compared to that of T . This delay suggests that a major damage to β-cells could occur in susceptible individuals before they test positive for islet-specific autoantibodies.
One could interpret panels (A1-C1) and (A2-C2) to correspond to high risk subjects that test positive to autoantibodies their entire lives but never develop the disease, panels (A3-C3) to correspond to high risk subjects that become type 1 diabetic and test positive to autoantibodies almost their entire life, and panels (A4-C4) to correspond to individuals that develop the disease very quickly (due to the presence of very potent and destructive T-cells), but test positive for autoantibodies during only a short window of time. Such criterion can be used to determine the risk associated with each clone of autoantibodies and the timing of T1D disease onset in individuals. S1.2. TWO-CLONE MODEL S1.2.1. Model Rescaling. We extend the one-clone model in Eqns. (S1.6)-(S1.11) by increasing the number of clones of T-cells under consideration. However, to limit the complexity of the model, we consider the simplest scenario possible in which we include two clones of T-cells with distinct autoantigenic specificities and different levels of avidity. Each one of these two clones is further divided into high and low avidity subclones. In other words, we assume that the model consists of two subclones T 11 and T 12 that are reactive to one autoantigen (labeled P 1 ) and two other subclones that are reactive to another autoantigen (labeled P 2 ), where k 22 ≤ k 21 ≤ k 12 ≤ k 11 .
Once again, we assume here that k 22 ≤ k 21 ≤ k 12 ≤ k 11 . The simulations presented in this paper are all generated using these equations and parameter values listed in Table S1.2.
(C4). The narrow regime with elevated level of β between the two black lines in (C3), representing the 0.3-critical threshold, is invoked by solution trajectories propagating close to the stable manifold of the transient saddle point. Although it is possible theoretically to target this regime to improve disease outcomes, its sensitivity to perturbations makes this clinically difficult, especially in later years when the regime is narrower.