Multiscale, mechanistic model of Rheumatoid Arthritis to enable decision making in late stage drug development

Model design and parameterization

The site of inflammation in RA is the joint, consisting of bone, cartilage, synovial membrane, and synovial fluid. Of these, the synovial membrane is the key site of inflammation with several immune cells and cytokines present¹⁵. We have assumed a well-mixed volume of 1 mL of synovium to be the site of interest to capture the basic pathophysiology of the disease. Inflammatory cells, cytokines and structural cells interact in our idealized synovial volume (Fig. 1 shows the modeling of Macrophages. Similar figures for all cell types and cytokines are available in Fig. 7). Once the model design process (more details in the Methods section) has resulted in the effect diagram and a set of equations for the model, the next step is determining the quantitative ranges of the parameters of the model. The broad approach is shown in Fig. 2, where “bottom-up” data from in vitro and preclinical experiments have been used to inform parameters and “top-down” data from clinical trials have been used to constrain the Vpop. We will focus on the Vpop in the Results and Discussion, but more details on model quantification have been discussed in the Methods section.

A single Calibrated Virtual Population consistent with available MTX, ADA and TCZ monotherapies was obtained

As described in the Methods section, we developed the model, quantified it, and generated a Vpop which is at dynamic ‘steady state’ with respect to the cells and cytokines represented in the model for various levels of severity of RA disease. Figure 3 shows the synovial cell densities of different cell types at baseline and in comparison, to reported ranges in the literature. (See Supplementary Data 1 for the literature from which these ranges were derived) The Vpop (n = 300) is well distributed within these ranges. We selected three therapies - Methotrexate (MTX), Adalimumab (ADA) and Tocilizumab (TCZ) to calibrate the Vpop from the Virtual Cohort, as discussed in the Methods. All patients in the Vpop were entered into the MTX and ADA trials and the subset of Virtual patients in this trial, which were non-responders to MTX (n = 251) were put on TCZ. Figure 3 shows the distribution of DAS28-CRP scores of the VPop at baseline, along with reported mean and standard deviation at baseline from the Kavanaugh et al. ¹⁶ and Yazici et al. ¹⁷ trials on ADA and TCZ respectively^16,17 (See Supplementary Data 1 for the detailed representation of available data from these trials). (The Strand et al. ¹⁸ study on MTX did not report baseline DAS28-CRP)¹⁸. The approach to derive disease scores from synovial densities is elaborated in the Methods section, under ‘Model parameterization: Implementation of disease severity scores’.

We implemented the PK and PD of the various interventions of interest - MTX, ADA and TCZ monotherapies as explained in the Methods section. In particular, the patients entering the MTX and ADA trials were “all comers” (i.e., not previously treated by either therapy but may have been treated by NSAIDS), whereas the TCZ trial was on MTX Inadequate Responders (IR). We simulated the trials as described starting with a single Vpop that is consistent with all the reported clinical scores reported (Fig. 4).

Figure 5 shows the results of calibration to therapy in terms of the metrics reported for clinical efficacy—ACR20, 50 and 70, as well as DAS28-CRP < 3.2 DAS28-CRP < 2.6. DAS28-CRP and dDAS28-CRP > 1.6. As can be seen from the figure, the Vpop calibration is a close match to data both for the MTX naive population and for the MTX-IR population used for calibration of TCZ. We also show that the distribution of DAS28-CRP scores in the VPop for each of these trials (see Fig. 11) shows sufficient phenotypic diversity in the response. The representation of DAS28-CRP scores in the model is discussed in the section “Additional information on methods, 4. Representing DAS28-CRP scores”. There can be tradeoffs in how to determine the Vpop fit, our goal was to have a range of poor and strong responders and thus directed our Vpop development accordingly.

An additional consideration in interpreting clinical scores in auto-immune trials is factoring in the response of the placebo arm in the trial. The efficacy attributable to the drug should be the placebo efficacy ‘subtracted’ from the reported efficacy. There are several approaches to implement placebo correction and we have followed the one recommended by Wang et al. ¹⁹, which is described in greater detail in the Methods section, under Stage 4 and 5¹⁹.

Validation: Predicting the outcome of a Tocilizumab trial on anti-TNF-α non-responder population

Model validation increases the confidence in the model – typically this is carried out by predicting impact of interventions which the model has not been calibrated to. This can be carried out by predicting the efficacy of novel therapies, combinations, or impact of calibrated therapies in new subpopulations. We validated our model by predicting the impact of Tocilizumab on a sub-population of MTX and ADA inadequate responders (IRs). This subpopulation (of n = 216) was selected as the set of Virtual patients with ACR < 50 and DAS28-CRP > 3.2 post-therapy for MTX and ADA. Figure 6 shows the outcome of our model simulations compared to the data obtained from a clinical trial on the efficacy of TCZ in such an ADA + MTX IR population of patients²⁰. As seen from Fig. 6, the predicted reduction in DAS28-CRP and ACR scores, with no further adjustments, using the calibrated Vpop shows a good match to ACR50, ACR70 and DAS28-CRP remission (i.e., DAS28-CRP < 2.6) scores. The match is poorer for ACR 20 though within 10%.

Stage 1: Establishing Project Goals

Our objective was to develop a model that could be used for the simulation of late-stage clinical trials (induction of Phase 2/3) that captures quantitatively the dynamics of cells and mediators in a mechanistic physiological scale and the impact of the therapy on the physiology and its translation to clinical disease activity scores within a time scale of a few months. We created a Vpop with varying disease severity of established, stable disease (vs. preclinical or early disease based on varying levels of cells and mediators in the joint)³¹. This Vpop is calibrated to be consistent with disease activity scores at baseline and in response to methotrexate and ADA, as reported in clinical trials^16–18,32. We aimed to develop a model structure that can add new biology and targets of interest within this framework; future versions will build on this scope and add functionality and be compatible with the features presented here.

Stage 2: Designing Model Scope and Complexity

The model design is as simple as necessary to capture the disease pathophysiology relevant to simulating the induction and maintenance phases of a typical Phase 2 or Phase 3 clinical trial, which runs for a few months in patients with established RA²⁹. The model is limited to capturing established steady-state disease (which can be mild, moderate, or severe) and does not capture the poorly understood complexities of disease onset, flares, and disease progression over several years^15,33.

The site of inflammation in RA is the joint, consisting of bone, cartilage, synovial membrane, and synovial fluid. Of these, the synovial membrane is the key site of inflammation with several immune cells and cytokines present¹⁵. We have assumed a well-mixed volume of 1 mL synovium to be the site of interest to capture the basic pathophysiology of the disease. Inflammatory cells, cytokines and structural cells interact in our idealized synovial volume.

The two types of cells considered in the synovial compartment of the model are structural and immune cells. In the synovial tissue, cell-cell contact-mediated effects, which may contribute to disease progression, are not considered at this time and all effects are assumed to occur via cytokines³⁴. The synovial volume itself is fixed over time in the model and the severity of disease is approximated by the cell density in the volume (higher for more severe disease).

Currently, other physiological components of interest, such as bone and cartilage are not captured in the model. Outside of the joint, we have a generalized plasma compartment. The lymph node is not modeled explicitly, and differentiated cells are assumed to migrate to synovial tissue from plasma. The model also does not consider heterogeneity between joints in a patient i.e., all joints of the patient are identical in their pathophysiology and response to therapy and only a single average joint is modeled. Cells and mediators involved in disease onset (such as Dendritic cells) and cells involved in other aspects of disease progression (like Osteoclasts involved in Bone erosion) were not captured in this version of the model. Autoantigen is not explicitly modeled, and we assume that excess autoantigen is always available (an argument for this approach is that tolerogenic therapies targeting autoantigen for RA have not shown much success³⁵.

A detailed description of the different cell types and cytokines in the synovial compartment of the model is provided below:

• A.

  Structural cells

• Fibrocytes like Synoviocytes (FLS): This is a specialized cell type located inside joints in the synovium and it is critical in the pathophysiology of RA disease. The inflamed synovium results in the formation of Pannus which is considered the characteristic feature of RA, and FLS are the most common cell type at the pannus–cartilage junction and contribute to joint destruction through their production of cytokines, chemokines, and matrix-degrading molecules and by invading joint cartilage^30,36.

• Endothelial cells (E): Endothelial cells play an important role in many diseases of chronic inflammation. In RA they express adhesion molecules and present chemokines, which leads to leukocyte migration from the blood into the tissue. Angiogenesis is another characteristic feature of RA and proliferating ECs lead to more blood flow and that leads to more trafficking of immune cells to the inflamed joint³⁷.

• B.

  Immune cells and Cytokines

In general, the inflamed synovium is characterized by infiltration of multiple immune cells and increased cytokine concentration in comparison with healthy tissue³⁸. At this time, we have included cells and mediators whose role in the disease is well-established and which have been considered for targeted therapies.

• Macrophages: Among the cells of the innate immune system, Macrophages are modeled. They are represented as a single homogeneous population, secreting both pro and anti-inflammatory cytokines. There is evidence that indicates greater complexity in the macrophages with multiple phenotypes in the inflamed RA joint, but this distinction may be made in later versions³⁸.

• Multiple Pro-Inflammatory Cells of the Adaptive Immune System: The cells of the adaptive immune system are central contributors to inflammation in RA³⁹. We have represented CD8 T Cells, Th1 and Th17 Cells, B Cells, and Plasma Cells.

• CD4 Treg: In RA, the role of anti-inflammatory cells is less well understood than that of pro-inflammatory cell types⁴⁰. However, Tregs are a key player that can induce strong mediator induced effects that may affect the dynamics of the other cell types⁴¹.

• Multiple Proinflammatory cytokines: Inflammation in RA is assumed to be primarily promoted by the pro-inflammatory cytokines which are numerous and often have redundant functionality^42,43. Several have been targeted in developing therapies with varying success. It is possible that the relative importance of cytokines may depend on the particular patient and their stage in disease progression³¹. We have included the following typical pro-inflammatory cytokines whose role in synovial inflammation has been studied and whose suppression has been successful in alleviating symptoms of stable RA disease in multiple studies: TNF-α, IL-6, IL-17, IL-12, and IL-23.

  Further, we have modeled IFN- γ and IL-1β as prototypical cytokines that have been shown to be important in establishing chronic inflammation. Some of their activity may be redundant with other cytokines, such as TNF-α, but their activity may also be necessary to incorporate as a ‘ceiling’ to the effectiveness of other anti-cytokine therapies. However, anti-IL-1β and anti-IFN-γ monotherapy have not been successful in mitigating symptoms of RA⁴⁴. GMCSF may affect disease via the activation of macrophages and has been modeled. Recent studies have shown some promise for anti-GMCSF therapy in mono or combination in alleviating symptoms of RA⁴⁵.

  B Cell Activating Factor, BAFF has been modeled to enable modulation effect on disease of B cell lineage⁴⁶. Auto-antibodies are also modeled and are considered to be biomarkers of inflammation in some patients as well to play an important role in cytokine secretion by immune cells⁴⁷. However, their impact on non-cytokine pathways (such as bone erosion, cartilage degradation) are not present. Thus, their effect on disease may not be captured comprehensively and can be improved in future versions.

  Not all cytokines with evidence of relevance to RA pathophysiology and treatment are included in the model. Cytokines such as IL-2 and IL-8, have not been explicitly included in the model; IL-2 is assumed to be present in non-limiting amounts in the synovium whereas IL-8’s main role is in the recruitment of PMNs (Poly Morphonuclear Neutrophils). PMNs are not captured in this version of the model since they are not thought to be major drivers of the disease and hence cytokines mainly relevant to PMN activity are also not explicitly captured³⁹. Cytokines such as IL-15, IL-18, IL-32 have shown some promise in improving clinical outcomes and may be very relevant to sub-populations of patients but have also not been included at this time. They can be added to the model structure and their functionality can be accounted for explicitly and quantitatively in subsequent versions. RANKL is another cytokine that is not included in the model that may have a major role to play in activating osteoclasts, which are not part of the current physiological scope.

• Multiple Anti-inflammatory cytokines: These cytokines down-regulate immune activation and reduce secretion of pro-inflammatory cytokines and are critical feedback to reduce inflammation but are inadequate in RA. We have modeled IL-10 and TGF-β, which are considered to be important anti-inflammatory mediators though, increasing these and other anti-inflammatory cytokines like IL-4 and IL-13 are generally insufficient to reduce disease symptoms in a population⁴⁰.

• C.

  Chemokines/ growth factors

  These are smaller proteins which are critical to immune cell migration to the joint and can also be important therapeutic targets⁴⁸.We have included a general factor regulating recruitment of all leukocytes that can be varied between virtual patients. This factor, called CAMS in the model (for Cell Adhesion Molecules) is at present a function of endothelial cells, to represent the expression of cell adhesion molecules on endothelial cells, and may be expanded further to denote specific adhesion factors⁴⁹. We have modeled RANTES (a key recruiter of T Cells), MIP-3α (recruits T Cells, B Cells and Macrophages), MCP-1 (recruits Macrophages). We also have modeled VEGF, a growth factor that induces the endothelial cells migration (which eventually increases angiogenesis) and exacerbates disease⁵⁰.

• D.

  Serum compartment

Elevated levels of cytokines are observed in both synovium and serum of RA patients. However, the measurements of cytokines in the serum tends to be highly variable and we have not attempted to match observations of serum cytokines to the model⁵¹. This central compartment in the model serves as a source of immune cells to be recruited into the synovial compartment. Further, the pharmacokinetics of therapy (MTX, ADA and TCZ) are modeled in this compartment and partitioned into the synovial compartment. Currently, we consider the serum compartment only as the input for the PK and as source of immune cells and have not attempted matching immune cell numbers in the serum, which are highly variable.

Stage 3: Developing Model Structure, Interaction Motifs among Model Components

Equations of Cells and Cytokines in the Synovial Compartment: The density of cells and concentration of cytokines is tracked in the synovium as a measure of disease severity (the synovium itself has a fixed volume of 1 mL in all Virtual Patients). Virtual Patients are modeled to have dynamic steady-state concentrations of cells and cytokines in disease which are reduced by therapy. The main motifs of interaction of the various constituents of the synovial compartment are: (i) Inflammatory cells release cytokines (ii) cytokines modulate life cycle of inflammatory and other cells (e.g., rates of proliferation, apoptosis and migration) and (iii) cytokines also modulate rate of secretion of cytokines by the inflammatory cells. The dynamics of each cell type’s density in the synovium is captured in an ordinary differential equation (ODE) and the following life-cycle processes are explicitly captured: proliferation, apoptosis, and migration into the joint from the central compartment (Fig. 1, and Fig. 7). Synovial cytokine concentration is captured by the model equations and tracks secretion of cytokines and other mediators by various cell types as well as their clearance. Each of the cell life-cycle processes are regulated by cytokines as shown in Fig. 1 (with macrophage life cycle as an example). We have assumed that the cells are well-matured, differentiated and activated in our model and these processes are not captured explicitly.

The dynamics of each cell is represented as shown in Eq. 1 below:

Where ProlifRate, MigRate and DegRate represent the rate of proliferation, migration or influx into the synovium and apoptosis rate of each cell type. Each of these processes is modeled as occurring at a baseline rate (cells/mL/hour) that can be further increased or decreased by cytokines or other mediators. Further, Proliferation and Migration rate are modeled as follows (Eqs. 2 and 3):

Where ProProlifEffect and AntiProlifEffect represent the fold increase and decrease in proliferation rate due to the effect of different cytokines that act on the cell. We have modeled the pro-inflammatory and anti-inflammatory effects as independent, with 0 < AntiProlifeffect < 1 and 0 < ProProlifEffect < 10. The effects of the cytokines on cell life-cycle rates are constrained based on available data, which is usually based on in-vitro experiments, and are saturating functions of cytokine concentrations across about 3 log scales. The equations describing these processes and design choices are explained in greater detail in the section ‘Additional information on Methods: 1. Representing cell and cytokine life cycle’.

Several approaches have been used to integrate cytokine signals which can be pleiotropic, redundant, synergistic, or antagonistic and may use overlapping intracellular pathways to affect the signaling^52,53. These include different ways of combining these signals in differential equation models such as ours as well as logic-based Boolean models^24,54–56. Our approach is designed to work at the scale of interest of reproducing the clinical effect of perturbations (anti-cytokine therapy) on a clinical population (vs., for example, reproducing cell cultures in an in vitro experiment)⁵⁷. More complex mechanisms such as increasing response of a cytokine by another cytokine (without a direct effect by the first cytokine on the process) or the involvement of obligatory cytokines for an²⁷ interaction are not modeled currently.

The dynamics of a cytokine in our model are captured in Eqs. 4–6 below:

Where Cytokine is the cytokine concentration, CytSecRate refers to the net cytokine secretion rate (ng/mL/hr) and CytDegRate (ng/mL/hr) refers to degradation rate of the cytokine.

CytDegRate is used from the literature where available, and generally ranges from minutes to hours. Data on cytokine half-life was obtained from multiple sources and often ranges from a few minutes to several hours (see Supplementary Data 2). The value BaselineRate_i refers to the baseline secretion rate by a particular cell type with density CellDensity_i. These rates are estimated to match reported steady-state concentrations of the cytokine. The relative amount secreted by each cell type was estimated to match data when available. ProSecEffect_i is the effect of the positive regulators on the baseline secretion rate and AntiSecEffect_i is the effect of the negative regulators of that particular cell type and are modeled similarly to their effect on cell life cycle rates. Greater detail on assumptions in the modeling are in the section, “Additional information on Methods: 2. Representing cytokine secretion and degradation”).

Stages 4 and 5: capturing behaviors, building confidence and exploring variability

Additional Information on methods

1. Representing Cell and cytokine Life cycle: The life cycle and regulation of the cells and cytokines represented in the model are shown in Fig. 7. The species in the model can be broadly classified into cells and soluble mediators. Each of these undergoes distinct processes. Cells undergo the following reactions: Proliferation, Influx and Degradation. The ODE for a cell type is written as follows:

Calculation of rates of proliferation, influx and degradation of the cell are described below.

1. Cell proliferation rate is calculated as:

   $$\mathit{ProlifRate}=\mathit{BaselineRate}*\left(1+\mathit{ProProlifEffect}\right)*(1-\mathit{AntiProlifEffect})$$

   9

   where

   ProlifRate is the rate of proliferation rate of each cell type, expressed in cells/(milliliter*day)

   BaselineRate is the baseline rate of proliferation, in absence of all modeled cytokines, expressed in cells/(milliliter*day)

   ProProlifEffect is the increase in the proliferation rate due to all cytokines which cause an increase in proliferation of that cell type.

   AntiProlifEffect is the reduction in the proliferation rate due to all cytokines which cause a decrease in proliferation of that cell type.

   We have chosen to model proliferation rate as a zeroth order term as it may be a more accurate representation of growth in a crowded, growth restricted environment, as well as allows the model to settle into a stable steady-state balancing the zeroth order proliferation and first order clearance. The combination of multiple cytokines is discussed in the next section and in the next section.

2. Cell migration or influx rate is calculated as:

   $$\mathit{InfRate}=\mathit{BaselineRate}*\left(1+\mathit{ProInfEffect}\right)*\left(1-\mathit{AntiInfEffect}\right)*(1-\mathit{AntiInfMTX})$$

   10

   where

   InfRate is the rate of influx rate of each cell type, expressed in cells/(milliliter*day)

   BaselineRate is the baseline rate of influx expected in absence of all modeled cytokines, expressed in cells/(milliliter*day)

   ProInfEffect is the proportional increase in the influx rate because of the cytokines/chemokines which enhance influx of each cell type.

   AntiInfEffect is the proportional reduction in the influx rate because of the cytokines/chemokines which reduce influx of each cell type.

3. Cell degradation rate is calculated as:

where

DegRate is the rate of removal of the cell type, expressed in cells/(milliliter*day)

BaselineRate is the baseline rate of degradation of the cell type expected in absence of all modeled cytokines, expressed in 1/day

CellDensity is the density of the cell type expressed in cells/milliliter

ProDegEffect is the proportional increase in the degradation rate of the cell type because of the cytokines which cause an increase in degradation of each cell type.

AntiDegEffect is the proportional reduction in the degradation rate of the cell type because of the cytokines which cause a decrease in degradation of each cell type.

• 2.
  Representing Cytokine secretion and degradation: Cytokines undergo the following reactions: Secretion and Clearance. The ODE for a cytokine is the sum of these reactions as follows:

  $$d\left(\mathit{Cytokine}\right)/\mathit{dt}=\mathit{CytSecRate}-\mathit{CytDegRate}$$

  12

  The secretion rate for each cytokine is a sum of secretion rates by several cell types and has the structure:

  $$\begin{array}{c}CytSecRate=\sum _i\left(BaselineRat{e}_i*CellDensit{y}_i*\left(1+ProSecEffec{t}_i\right)\right)\hfill \\ *\left(\left(1-AntiSecEffec{t}_i\right)\right)\hfill \end{array}$$

  13

  where

  CytSecRate is the total rate of secretion of the cytokine, expressed in nanogram/(milliliter*day)

  BaselineRate is the baseline rate of secretion of the cytokine by a cell type expected in absence of all modeled cytokines, expressed in nanogram/(cell*day)

  CellDensity is the density of each cell type expressed in cell/milliliter

  ProSecEffect is the proportional increase in the secretion rate of the cytokine by a cell type because of the effect of other cytokines on the cell.

  AntiSecEffect is the proportional reduction in the secretion rate of the cytokine by a cell type because of the effect of other cytokines on the cell.

  The clearance rate of a cytokine is expressed as follows:

  $$\mathit{CytDegRate}=\mathit{BaselineRate}*\mathit{CytokineConc}$$

  14

  where

  CytDegRate is the rate of removal of the cytokine, expressed in nanogram/(milliliter*day)

  BaselineRate is the baseline rate of clearance of the cytokine expressed in 1/day

  CytokineConc is the concentration of the cytokine expressed in nanogram/milliliter

• 3.
  Representing effect of multiple cytokines: As observed above, the rates of reactions have been moderated using terms such as:

  $$\left(1+ProRate\right)\mathrm{and}\left(1-AntiRate\right)$$

  These refer to the increase and decrease in the cell life-cycle rates because of the complex cytokine milieu of the synovium. Each term is a sum of individual contributions from several cytokines and has the following structure:

  $$\mathit{ProRate}=\min (\mathit{Prolimit},\sum _{\mathit{cytokines}}\mathit{MM}(\mathit{CytokineConc},\mathit{VmEffect},\mathit{KmEffect},\mathit{SlopeEffect}))$$

  15

  $$\mathit{AntiRate}=\min (\mathit{Antilimit},\sum _{\mathit{cytokines}}\mathit{MM}(\mathit{CytokineConc},\mathit{VmEffect},\mathit{KmEffect},\mathit{SlopeEffect}))$$

  16

  where

  ProRate is the proportional increase in any rate because of the surrounding cytokines

  Prolimit is the limit imposed on the proportional increase of any rate to accommodate biological feasibility; currently, this is fixed at 10-fold.

  Antilimit is the limit imposed on the proportional decrease of any rate to accommodate biological feasibility; currently this is fixed at 0.75.

  MM is an external function that calculates the value of the Hill expression

  $$\frac{\mathit{VmEffect}*{\mathit{CytokineConc}}^{\mathit{SlopeEffect}}}{{\mathit{CytokineConc}}^{\mathit{SlopeEffect}}+{\mathit{KmEffect}}^{\mathit{SlopeEffect}}}$$

  VmEffect is the maximum possible increase in a rate by the cytokine

  KmEffect is the concentration of the cytokine concentration at which the rate is increased by half the maximum possible increase by the cytokine

  SlopeEffect is the hill coefficient of the hill function

  CytokineConc is the concentration of the cytokine expressed in nanogram/milliliter

  Both Pro and Anti effects are capped to lead to a maximum 10-fold change in the activity. But they hold different values by nature of how they are used in the equations. For example, ProProlifeffect is used as “Rate = Baseline*(1+ProProlifeffect)” and a cap of 10 on ProProlifEffect leads to a maximum of “Rate = 11*Baseline”. Meanwhile, AntiProlifeffect is used as “Rate = Baseline(1-AntiProlifeffect)” and a cap of 0.9 on AntiProlifeffect leads to a minimum of “Rate = 0.1*Baseline”. Both lead to a similar fold change from baseline, but in opposite directions.

• 4.
  Representing DAS28-CRP score: The following expression is used to calculate the DAS28-CRP value of the virtual patients.

  $$\mathit{DAS}28-\mathit{CR}{P}_{\mathit{model\; score}}=2.5*\mathit{Hill}\left(\mathit{Macrophages}\right)+2*\mathit{Hill}\left(\mathit{FLS}\right)+1.5*\mathit{Hill}\left(\mathit{Th}1\right)+1.5*\mathit{Hill}\left(\mathit{BCell}\right)+1*\mathit{Hill}\left(\mathit{PlasmaCell}\right)+0.5*\mathit{Hill}\left(\mathit{Th}17\right)+0.5*\mathit{Hill}\left(\mathit{Endothelial}\right)+0.5*\mathit{Hill}\left(\mathit{CTL}\right)-0.5*\left(\mathit{Treg}\right)$$

  17

  Where

  $$\mathit{Hill}\left(\mathit{Cell}\mathrm{\_}x\right)=\left(\mathit{Cell}{\mathit{Density}\mathrm{\_}x}^{\mathrm{\wedge }}\mathit{gamma}\mathrm{\_}x\right)/({\mathit{Cell}\mathrm{\_}\mathit{Density}\mathrm{\_}x}^{\mathrm{\wedge }}\mathit{gamma}\mathrm{\_}x+{\mathit{Km}\mathrm{\_}x}^{\mathrm{\wedge }}\mathit{gamma}\mathrm{\_}x)$$

  The values of Km and gamma for each cell type are in Table 2 below.

  The Km_x for each cell_x is calibrated such there is a uniform distribution of $\mathit{Hill}\left(\mathit{Cell\_x}\right)$ from 0 to 1, at baseline for the cohort and that eventually, overall DAS28-CRP is well distributed between 0 and 10. The gamma_x was also chosen to be at 2.5, to achieve the mentioned behavior. The gamma_x is fixed to be the same across all cell types such that similar fold reductions in the cell densities translate to similar reductions in their subscores. Once the DAS28-CRP score has been calibrated, Virtual Patients in the cohort with baseline DAS28-CRP < = 3.2 were removed from the cohort as they correspond to low disease activity and are not eligible for the clinical trials.

• 5.
  Sensitivity Analysis of the Model: We carried out both local and global sensitivity analyses on the parameters of the model to understand the dependence of the model outcome on the different parameters of the model.

  1. Local Sensitivity Analysis: From the VPop, the virtual patient closest to the median synovial cell densities was picked for this analysis. To measure the sensitivity of each parameter, the percent change in DAS28-CRP score is observed by introducing 2x variability to the parameter value of the virtual patient. The top 20 sensitive parameters are plotted as a tornado plot as shown in Fig. 8. As can be seen from the figure, macrophage related parameters are the most sensitive ones in this patient, along with parameters affecting influx of all immune cells (Leukoinflux_MaxbyCAM), FLS proliferation (kg_FLS_baseline), Th1 proliferation etc. This is in line with our expectations that macrophages are important drivers of the disease. Note that these sensitivities are likely to be different in other patients in the Vpop; in fact, this is desirable since we want a diversity of responses to therapy, driven by a diversity of underlying mechanisms of disease. The limitation of the local sensitivity analysis is that the sensitivity observed to the parameters may be very specific to the virtual patient/parametrization around which the variability is generated. This means that the local sensitivity analysis may not represent the sensitivity of the model across the whole parameter space explored.
  2. Global Sensitivity Analysis: To address the shortcomings of the local sensitivity analysis, we have also run a global sensitivity analysis. The first order and total order Sobol indices have been calculated for each parameter using a sample size of 10000 within the explored parameter space. The total order indicates the contribution of a particular parameter to total variance in the outcome; including its impact due to joint parameter variations while the first order sensitivity indicates sensitivity to the particular parameter alone⁹⁶. The top 20 sensitive parameters have been plotted in Fig. 9. As can be seen from the figure, FLS, macrophage and B cell related parameters are the ones which show the most impact on disease severity; the parameter F-CAM regulates the influx of all cell types in the model and is a lumped representation of multiple chemokines. Total order is greater than first order sensitivity for most of the parameters, indicating that co-variance/interactions between the parameters contributes more to model variability than variability in each single parameter.
• 6.

  Pharmacokinetics and Pharmacodynamics of therapies: The pharmacokinetics of the drugs used are modeled and parameterized based on the existing work in public literature. The PK models of Methotrexate, Adalimumab and Tocilizumab are implemented as described as summarized in Table 3.

1. Methotrexate PD: Methotrexate is modeled to affect the system in the model, in the following three routes^89–92.

1. Reduce the secretion rate of cytokines by all cells except Tregs.

2. Increase the production rate of Anti-Inflammatory cytokines by Tregs.

3. Reduce the influx rate of Immune cells into the synovium.

   These effects are calculated using the following Hill expressions, respectively:

   $$\begin{array}{c}Ant{i}_{CytSe{c}_{MTX}}=MM\left(MTXDrug\_Central\_available,Anti\_CytSec\_MaxbyMTX,\right)\hfill \\ \left(HalfEffectConc\_Anti\_CytSec\_byMTX,Slope\_Anti\_CytSec\_byMTX\right)\hfill \end{array}$$

   18

   $$\begin{array}{c}Pro\_CytSec\_MTX=MM\left(MTXDrug\_Central\_available,Pro\_CytSec\_MaxbyMTX,\right)\hfill \\ \left(HalfEffectConc\_Pro\_CytSec\_byMTX,Slope\_Pro\_CytSec\_byMTX\right)\hfill \end{array}$$

   19

   $$\begin{array}{c}Pro\_CytSec\_MTX=MM\left(MTXDrug\_Central\_available,Anti\_CellInflux\_MaxbyMTX,\right)\hfill \\ \left(HalfEffectConc\_Anti\_CellInflux\_byMTX,Slope\_Anti\_CellInflux\_byMTX\right)\hfill \end{array}$$

   20

   Where
   $\mathit{MM}(\mathit{DrugConc},\mathit{VmEffect},\mathit{KmEffect},\mathit{SlopeEffect})$ is an external function that calculates the value of the Hill expression

   $$\frac{\mathit{VmEffect}*{\mathit{DrugConc}}^{\mathit{SlopeEffect}}}{{\mathit{DrugConc}}^{\mathit{SlopeEffect}}+{\mathit{KmEffect}}^{\mathit{SlopeEffect}}}$$

   $\mathit{VmEffect}$ is the maximum possible proportional change in a rate by the drug.

   $\mathit{KmEffect}$ is the concentration of the cytokine concentration at which the rate is changed by half the maximum possible change by the drug.

   $\mathit{SlopeEffect}$ is the hill coefficient of the hill function.

   Vm, Km and SlopeEffect are calibrated at 0.5, 1E-5 mg/L and 2, respectively, to achieve a reasonable mean reduction in disease severity score with therapy.

   $\mathit{MTXDrug\_Central\_available}$, expressed in milligram/liter, is the central concentration of methotrexate multiplied by the bioavailability of the drug.

   $\mathit{Anti\_CytSec\_MTX}$ is the proportional reduction in the secretion of the proinflammatory cytokines

   $\mathit{Pro\_CytSec\_MTX}$ is the proportional increase in the secretion of the anti-inflammatory cytokines

   $\mathit{Anti\_CellInflux\_MTX}$ is the proportional reduction in the influx of the immune cells.

   The therapy effect is achieved by multiplying the corresponding cytokine secretion rates by pro-inflammatory cells, cytokine secretion rates by Tregs and the Influx rates of cells into synovium by $(1-\mathit{Anti\_CytSec\_MTX})$, (1 + $\mathit{Pro\_CytSec\_MTX})$ and $(1-\mathit{Anti\_CellInflux\_MTX})$, respectively.

   1. Adalimumab PD: The effect of adalimumab is modeled as a reduction in the free TNF-α, due to the formation of an ADA_TNF- α complex, after binding to the drug in synovium.
   2. Tocilizumab PD: Tocilizumab targets the IL-6 receptor. As there is no explicit representation of IL-6 receptor in the model, we have assumed that a similar PD effect can be achieved by reducing the IL-6 by a similar proportion as the receptor is expected to be. This is implemented by dividing the clearance of IL-6 with $KD\_TCZ/\left(TCZDrug\_Synovium+KD\_TCZ\right).$

The above factor corresponds to the fold reduction expected in IL-6 receptors, in a scenario where the Drug concentration is significantly higher than the receptor IL-6 concentration.

• 7.
  Development and characterization of the Vpop:

  1. The Vpop was generated using a standardized workflow depicted in Fig. 10.
  2. The distribution of clinical scores before and after therapy was determined as shown in Fig. 11to confirm that the Vpop is not biased and has a broad distribution of clinical scores.
  3. Variability was introduced to a subset of the model parameters to generate the virtual population. Figure 12 describes the distributions of these parameters in the virtual population on a logarithmic X-axis, while Table 4 shows the ranges of the parameters shown in Fig 12 and provides a key to read Fig 12.
  4. The distribution of cells and cytokines in the model at steady state was also determined to show that the entire range of cell and cytokine values obtained from the literature were spanned by the Vpop, as seen in Figs. 13 and 14. Table 5 provides the log transformed ranges of the cell densities at baseline in the Vpop while Table 6 provides the log transformed ranges of cytokine concentration at baseline in the Vpop.

Software details

The model consists of a series of ODEs representing the dynamics of the cells, cytokines and regulatory interactions represented. The general form of these ODEs is described in the section on Model reactions and in the section, “Additional information on Methods”. The model was created and simulated in MATLAB Simbiology. The virtual cohort generation and simulations were run using gQSPsim⁹⁴. The Vpop selection was done using a genetic algorithm developed by Vantage Research. Plots were generated using MATLAB.