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

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,17,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:

1. 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³⁷.

1. 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⁴⁰.

1. 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⁵⁰.

2. 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 ProProlifEffect

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,55,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”).

Model Parameterization

In all cases, data is taken from “bottom-up” parameters which determine the rates of the various reactions and is used to directly inform the model parameters. “Top-down” data from clinical, synovial, RA samples is prioritized over clinical sera or pre-clinical measurements to determine model behavior and further constrain the parameter space. Below, we discuss in further detail the sources of parameterization for each subset of model parameters. Documentation of all model parameters estimated from literature and used in the model is available in the Supplementary Data 2.

When possible, ‘bottom-up’ data from in vitro experiments are used to determine the rate constants associated with cellular processes such as apoptosis, migration etc. Mediator driven effects are also determined from in vitro experiments, but the model quantification is limited to ‘reasonable’ fold changes in vivo. For example, the fold increase in proliferation due to cytokines is capped at 10x of baseline rates. When such parameters are unavailable, they are estimated such that the model is consistent with ‘top-down’ data such as measured synovial cell densities and cytokine concentrations at steady-state and in response to therapy (if such measurements are available on therapy). Steady state cell numbers for all cell types used in the model are obtained from literature data on cell numbers in terms of cells/mm² from IHC/ microscopy based examination of biopsy samples of RA patients and converted them into cells/mL based on the method described in Preza et al. ⁵⁸ The experiments and related calculations are described in the “Supplementary Data 1” in the sheet “conversion_ cells_ml”.

In particular, the following are some examples of rate constants of cell life cycle and cytokine secretion used to parameterize the model:

• Rate of cell apoptosis (1/hour): This was calculated from the literature data on the % apoptosis seen in activated immune cells (primary cell lines like T Cells, Macrophages etc. derived from healthy human PBMC, usually treated with LPS (vs. passaged cell lines)(See for e.g. ⁵⁹). Among structural cells, for the endothelial cells, half-life measured in inflamed tissue was used⁶⁰. For FLS, the % apoptosis of FLS reported in in-vitro cell cultures of RA -FLS was used⁶¹. Apoptosis rates for various cell types in in-vitro settings are reported in the literature and tend to range from a few hours to a few days. For e.g., to calculate the half-life of Th1-cells, we used data on activation induced cell death (AICD) i.e., % apoptosis reported in activated Th1 cells when stimulated with anti-CD3⁶². This data showed 15% apoptosis in 6 hours which is then translated into a half-life of 0.65/day using the equation t_1/2= ln (100/ (100-% apoptosis)/time of measurement in days.

• Rate of cell migration or influx into synovium (cell/ml/hour): These rates, when available were determined from migration assays carried out for the different cell types e.g., for determining the baseline (i.e., unregulated) migration rate of CD4 T Cells into the synovium, we looked at an in-vitro migration assay where the fraction of CD4 T cells that cross a semi- permeable membrane in absence of any chemokine in 2 hours was noted⁶³. This rate is then increased by different folds in response to different chemokines.

Two different processes contribute to increase in cell density of any cell type in the model – migration or influx of the cell type into the synovium in response to chemokines and proliferation in response to different signals received by the cell type. The relative contribution of synovial cell density due to proliferation within the synovium vs. influx from the serum is poorly characterized in the literature. For e.g., a study tracked autologous labeled macrophages in RA patients and found 64. No such data was available for lymphocytes. In steady-state disease, in the Vpop, the flux due to proliferation tends to dominate across cell types and Virtual Patients in the model and is consistent with our best understanding based on clinical trial data where anti chemokine therapies have not been successful in RA whereas anti proliferative therapy e.g., Abatacept (Anti-CTLA4) has been successful^65,66.

• Rate of proliferation (cell/ml/hour): Given reasonable constraints on the other 2 fluxes of cell life cycle, rates of proliferation were determined by fitting to obtain the observed steady state number of the cell type in RA. The proliferation rate is assumed to be of zeroth order, as is needed to maintain a steady state of the disease.

• Rate of cytokine secretion (ng/cell/hour): For calculating cytokine secretion rate per cell type in the reference Virtual Patient, we have used the following strategy:

• The concentration of cytokine in the tissue (synovial tissue or synovial fluid) is determined from the literature (see the Supplementary Data 2 for references).

• The relative secretion of cytokine by the various cells is determined from literature. For e.g., TNF-α can be secreted by macrophages, Th1 cells, CD8 T cells, B cells and FLS in the model. From literature, we determine that macrophages are the major contributor of the TNF-α seen in the joint, while T-cells are also known to be important secretors of TNF-α^67,68.

• Cytokine degradation rate is obtained from the reported half-life of each cytokine (see Supplementary Data 2). The cytokine secretion rates were varied during the creation of the virtual cohort and the relative order of secretion rates is not maintained in the final Vpop.

In addition, rates of each of these processes can be regulated by cytokines. Regulation of cell behavior by a cytokine is assumed to be a fold-increase over a baseline rate and is a saturating function of cytokine concentration. In some cases, a particular cytokine may be essential for a behavior e.g., IL-12 is essential for IFN-γ secretion by Th1 cells⁶⁹. In this case, we do not represent the effect as a fold increase, but rather modify the function to make IFN-γ secretion 0 in the absence of IL-12. The maximum effect of a cytokine in increasing or decreasing a rate is capped using a value determined from the literature (see Supplementary Data 2 for the quantification of these values). In the case of multiple cytokines regulating the same effect, the total effect is capped at 10-fold increase or 10-fold decrease from baseline rates (See the section “Additional information on Methods, 3, Representing the effect of multiple cytokines” for the mathematical representation of this). This is a heuristic which we hypothesize avoids unphysiological behavior clinically, though in vitro data may show higher increases. The variability in patient response may be sensitive to this hypothesis and must be investigated further.

Implementation of disease severity scores: Once the basic interactions of the cells and cytokines in the synovium are set up, this physiological scale is connected to clinical disease status. Multiple disease severity scores, such as the DAS28-CRP, ACR, SDAI have been tracked in clinical trials as markers of baseline disease characteristics and response to therapy^70,71. These measures track subjective and objective assessment of disease severity in RA and do not require measurement of immune cells or mediators from serum or synovium. Our assumption is that the densities of cells driving the disease correlate with disease severity and hence the disease score. This step of connecting the cell densities in the model to disease severity scores, which are based on subjective assessments, while reasonable, is an empirical assumption attached to an otherwise mechanistic model.

Other models have used similar approaches²⁴. This connection between disease score and actual disease severity as measured by serum and synovial measurements, imaging of affected joints have been explored with other biomarker studies (Correlation of disease score with RF, anti CCP, serum cytokine levels^71,72. In particular, a multi-biomarker panel based on serum measurements was developed and validated to (VECTRA^73,74,75), other serum proteins. The connection between biomarkers and clinical scores in RA needs further investigation to establish a stronger connection between physiology and disease severity^76,77,78.

The model is calibrated to DAS28-CRP and ACR as the clinical scores. DAS28-CRP has been shown to have robust correlations with disease and eventual patient outcomes of quality of life and is being reported in most recent trials^79,80,81. Quantitatively, one major advantage is that DAS28-CRP scores are reported at baseline as an entry criterion (vs. the ACR score which is only defined as a change in a study). In developing the expression connecting the synovial read-outs to DAS28-CRP score, using our best understanding of the disease physiology and prominence of different cells to the disease activity, we have allocated different coefficients to different cell types to contribute to the final score. The other parameters such as Kms and γs were calibrated based on the physiological steady state ranges for their corresponding cell types.

Where

The values of Km and γs for different cell types are provided in Table 2 along with the methodology using which they were calibrated (See ‘Additional information on Methods: 4. Representing DAS28-CRP). As mentioned earlier, the ACR score does not have a baseline and is tracked as a change from baseline. For this reason, we have calculated the ACR score in our model as a percent change in the DAS28-CRP, from baseline. We made this as the simplest assumption given limited understanding of the physiological differences between these scores.

The DAS28-CRP in the model is driven by the best available constraints and modeler choices to maintain a simple, transparent structure. Other modelers may be motivated to derive a more intricate mathematical model to connect synovial species and clinical scores, and this can be a major area of model refinement if more data is available in proprietary or other data sources.

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

Characterizing steady-state behavior

Cell numbers were derived from synovial biopsies of RA patients with established disease of moderate to severe severity while cytokine values were derived either from synovial biopsies or synovial fluid samples. We derive the ranges for each cell type/cytokine from the data reported in either synovium or synovial fluid (see Supplementary Data 1). Further we have performed Local and Global Sensitivity Analyses and have identified key model parameters that influence the disease state as shown in Figs. 8 and 9 (Also see ‘Additional information on methods: 5. Sensitivity Analysis of the model). We ensured that all the virtual patients used in Vpop calibration are constrained to lie within the range derived from literature.

The figure shows the percent change in DAS28-CRP score observed when introducing 2x variability to the parameter values of the virtual patient with the blue bars representing the outcomes of a 2x increase and the red bars representing the outcomes of a 2x decrease in the value of the parameter. E.g., on increasing the parameter MacroApop_MaxbyIFng from its baseline value by 2x, there is a 25% increase in the DAS28-CRP score while a 2x reduction in the same parameter leads to an 8% reduction in the DAS28-CRP score.

The total order (represented by the blue bars) 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 (represented by the red bars) indicates sensitivity to the particular parameter alone. E.g., we can see that kg_FLS_Baseline contributes maximally to the total variation seen in the DAS-28CRP score, both directly and through its interactions with other parameters.

A reference virtual patient (ref VP) was generated by calibrating the following parameters to replicate the median cell and cytokine behavior.

1. 1.

   Proliferation of Cells

2. 2.

   Secretion rate of cytokines

Once the reference virtual patient was calibrated, we generated multiple virtual patients, i.e., a virtual cohort, (n ~ 200000) by varying the following parameters.

1. 1.

   Proliferation and migration rates of cells

2. 2.

   Secretion rate of cytokines

3. 3.

   Upregulation/Downregulation of cytokines on cell processes

The cell ranges from literature are used as constraints on the generated virtual patients to determine which patients are ‘physiologically plausible’ and will be included during calibration to therapies. This set of plausible virtual patients (n = ~50000) is termed as a Virtual Cohort. While developing this virtual cohort, we have given greater importance to the values of cells than to the values of cytokines determined from the literature. This is because cytokine measurements are generally not from the tissue itself (some are from synovial fluid, while others are from plasma) whereas cell measurements are from the inflamed tissue. It is at this stage that we additionally constrain the virtual patients such that their DAS28-CRP is >3.2.

Characterizing response to therapy

After virtual cohort generation to ensure reasonable behavior of virtual patients at steady-state i.e., in absence of therapy, we calibrated the response of the virtual patients to different therapies via Vpop generation. The different criteria based on which we selected therapies for calibration are discussed in the section below (Selection criteria for trials for model calibration). In some cases, a particular cytokine may be essential for a behavior e.g., IL-12 is essential for IFN-γ secretion by Th1 cells⁶⁹. In this case, we do not represent the effect as a fold increase, but rather modify the function to make IFN-γ secretion 0 in the absence of IL-12. Based on these criteria, we selected Methotrexate (MTX), Anti-TNF-α and Anti-IL-6R as the drugs to which the response of the Vpop is matched. We selected MTX therapy to model since this is considered as an ‘anchor drug’ for RA^79,80 and most patients entering clinical trials for other therapies display an inadequate response to MTX or continue to be on MTX therapy to control RA symptoms^82,83. When patients show an inadequate response to MTX, biologic therapy is often the next line of treatment. Anti-TNF-α is the first biologic prescribed to RA patients who have inadequate response to DMARDS such as MTX. Five Anti-TNF-α have been approved for RA so far – Infliximab, Etanercept, Adalimumab, Certolizumab Pegol and Golimumab⁸⁴. All of them have similar efficacy in terms of remission and response but have different modes of delivery (Infliximab is an IV drug while the others are given SC) and different half-lives⁸⁵. Of these, Adalimumab was selected for model calibration. Tocilizumab or Anti IL-6R was included into the model calibration since it is one of the therapies (along with Rituximab or Abatacept) that are tried out in anti-TNF-α IR patients⁸⁴. The model also has adequate structure to implement novel therapies like Anti-IL-17, Anti-IL-1 and JAK/Tyk inhibitors etc. that are developed to address non-responders to the first line therapies.

Selection criteria for trials for model calibration: The design and purpose of the model needs to be considered carefully in selecting the trials to calibrate the Vpop against. The clinical trials for the therapies to be calibrated are selected using the following criteria:

• Steady-state, induction phase trials: In this case, the model developed is of steady-state disease (for 12 months) and will need to include disease progression which is not a feature of this model.

• Trial design: Trials with placebo arms are prioritized; in absence of any trials with a placebo arm, a trial with a common comparator of already existing therapy can be used. Trial protocols should be carefully examined for any changes such as crossover design which can complicate the interpretation of results. Phase III trials are generally preferred because of the large size of the sample.

• Patient population: Different population characteristics are of interest to the community and should be considered when calibrating the Vpop. Anti-TNF-α therapy is the recommended first biologic therapy for inadequate responders to DMARDS. Hence, anti-TNF-α non-responders and anti-TNF-α naive patients are explicitly recruited, in RA clinical trials for novel therapies to determine the response of both types of patients to novel therapies (e.g., see¹⁷. The Vpop generated from the model should contain both these sub-populations.

• Trial outcomes: It should be noted that clinical trials of interest should have the clinical readout of interest. For example, if ACR and DAS28-CRP are not reported in a clinical trial paper, it may not be preferred. For example, in the context of RA, the distribution of DAS28-CRP scores at baseline, change in DAS28-CRP mean value at different time points at different doses, remission and response scores should be recorded and investigated for conflict or concord.

Treatment of Placebo: The gold standard for judging the efficacy of a novel therapeutic is the randomized, controlled, clinical trial where the remission and response seen in the treatment arm should be compared to that seen in the control arm⁸⁶. Patients with autoimmune disease can show spontaneous response/remission, hence, a quantitative comparison of therapy with placebo as the comparator can be useful to determine efficacy. However, placebo correction by absolute subtraction of placebo arm remission from treatment arm remission does not account for the possibility that patients who respond to placebo can also respond to therapy. This can lead to an underestimation of the true therapeutic efficacy of the drug. To account for this, Wang et al. ¹⁹ uses a method of probabilities to arrive at a placebo correction method, which we have used in our model¹⁹.

The summary of trial characteristics of selected trials for model calibration are shown in Table 1.

It is to be noted that the trials for Adalimumab and Tocilizumab were run in combination with Methotrexate and with Placebo + Methotrexate as the control arm. Due to this, we have considered the placebo corrected clinical outcomes to represent the efficacy of Adalimumab and Tocilizumab monotherapies and used for the calibration of monotherapies in the model.

Modeling PK and PD of the therapies

The model captures the PK of the drug in the tissue and serum levels using the simple two or three compartment models. The PD of the drug is captured by using a binding kinetics model relevant to the drug.

1. 1.

   PK modeling: The parameters related to PK of the therapies in the model were taken from literature and are described in the section, “Additional information on Methods, 6. Pharmacokinetics and Pharmacodynamics of therapies”.

2. 2.

   PD modeling:

1. a.

   Anti-cytokine therapies such as Adalimumab are modeled using the in-vitro binding parameters assuming equilibrium kinetics^87,88.

2. b.

   Anti-Cytokine Receptor therapy such as Tocilizumab is modeled with an assumption that the reduction in the free receptor has equivalent effect as similar fold-reduction in free cytokine concentration. This assumption is made because the cytokine functions by binding to its receptor and hence blocking the receptor should have equivalent effects as blocking the cytokine itself.

3. c.

   MTX therapy affects (a) cytokine secretion (b) migration and (c) Treg function^89,90,91,92. The effect of the dose-response in each pathway is adjusted to match the clinical readout. The equations representing the PD are in the section, “Additional information on Methods, 6. Pharmacokinetics and Pharmacodynamics of therapies”.

Development of Virtual Population

The selected patients were then subjected to the therapies to which the Vpop is to be calibrated (here, MTX, Ada and TCZ). A subset of the virtual cohort, which shows appropriate baseline characteristics and responses to therapies – as determined by comparison with the values noted from the clinical trials, is selected to form the final Vpop. A Vpop of ~300 VPs were selected using this method. This process is described in Fig. 10. The distribution plots for the model parameters in the Vpop and the distribution plots of resultant steady state cell and cytokine concentrations are included in the section, “Additional information on Methods, 7. Development and characterization of Vpop”.

Flow chart depicting the process of creation of a virtual cohort and its calibration to multiple therapies to generate a Virtual population.

Cohort Enrichment: We had to enrich a few patient phenotypes in the cohort before we were able to generate a reasonable Vpop. Enrichment involves recognizing the under-represented phenotypes in the cohort that need to be more prevalent to obtain a good Vpop and generating similar plausible virtual patients by re-sampling from a tighter parameter space around them⁹³. This process can be iterative. To make the cohort enrichment efficient, we used prevalence weighted Vpop generation using gQSPsim and used the patients with very high prevalence weights to be enriched⁹⁴. Other groups have developed QSP models with multiple Vpops in order to improve the coverage of the plausible parameter space and offer even more robust predictions^26,95. Developing multiple VPops for this model is definitely a future area of improvement.

Additional Information on methods

1. 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. (a)

   Cell proliferation rate is calculated as:

   $${ProlifRate}={BaselineRate}* \left(1+{ProProlifEffect}\right)* (1-{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. (b)

   Cell migration or influx rate is calculated as:

   $${InfRate}={BaselineRate}* \left(1+{ProInfEffect}\right)* \left(1-{AntiInfEffect}\right)* (1-{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. (c)

   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.

1. 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({Cytokine})/{dt}={CytSecRate}-{CytDegRate}$$

   (12)

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

   $$\begin{array}{l}{CytSecRate}=\mathop{\sum }\limits_{i}\left({BaselineRat}{e}_{i}* {CellDensit}{y}_{i}* (1+{ProSecEffec}{t}_{i})\right.\\\qquad\qquad\qquad\;\;* \left.(1-{AntiSecEffec}{t}_{i})\right)\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:

   $${CytDegRate}={BaselineRate}* {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

2. 3.

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

   $$\left(1+{ProRate}\right)\,{\rm{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:

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

   (15)

   $${AntiRate}=\min ({Antilimit},\sum _{{cytokines}}{MM}({CytokineConc},{VmEffect},{KmEffect},{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{{VmEffect}* {{CytokineConc}}^{{SlopeEffect}}}{{{CytokineConc}}^{{SlopeEffect}}+{{KmEffect}}^{{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.

3. 4.

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

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

   (17)

   Where

   $${Hill}({Cell}{\rm{\_}}x)=({Cell}{{Density}{\rm{\_}}x}^{{\rm{\wedge }}}{gamma}{\rm{\_}}x)/({{Cell}{\rm{\_}}{Density}{\rm{\_}}x}^{{\rm{\wedge }}}{gamma}{\rm{\_}}x+{{Km}{\rm{\_}}x}^{{\rm{\wedge }}}{gamma}{\rm{\_}}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 \({Hill}({Cell\_x})\) 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 

4. 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. a.

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

      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.

5. 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. a.

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

1. 1.

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

2. 2.

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

3. 3.

   Reduce the influx rate of Immune cells into the synovium.

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

   $$\begin{array}{ll}{Ant}{i}_{{CytSe}{c}_{{MTX}}}={MM}\left({MTXDrug}{{\_}}{Central}{{\_}}{available},{Anti}{{\_}}{CytSec}{{\_}}{MaxbyMTX},\right.\\\left.\qquad\qquad\qquad\quad{HalfEffectConc}{{\_}}{Anti}{{\_}}{CytSec}{{\_}}{byMTX},{Slope}{{\_}}{Anti}{{\_}}{CytSec}{{\_}}{byMTX}\right)\end{array}$$

   (18)

   $$\begin{array}{ll}{Pro}{{\_}}{CytSec}{{\_}}{MTX}={MM}\left({MTXDrug}{{\_}}{Central}{{\_}}{available},{Pro}{{\_}}{CytSec}{{\_}}{MaxbyMTX},\right.\\\left.\qquad\qquad\qquad\qquad\qquad{HalfEffectConc}{{\_}}{Pro}{{\_}}{CytSec}{{\_}}{byMTX},{Slope}{{\_}}{Pro}{{\_}}{CytSec}{{\_}}{byMTX}\right)\end{array}$$

   (19)

   $$\begin{array}{ll}{Pro}{{\_}}{CytSec}{{\_}}{MTX}={MM}\left({MTXDrug}{{\_}}{Central}{{\_}}{available},{Anti}{{\_}}{CellInflux}{{\_}}{MaxbyMTX},\right.\\\left.\qquad\qquad\qquad\qquad\qquad{HalfEffectConc}{{\_}}{Anti}{{\_}}{CellInflux}{{\_}}{byMTX},{Slope}{{\_}}{Anti}{{\_}}{CellInflux}{{\_}}{byMTX}\right)\end{array}$$

   (20)

   Where

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

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

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

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

   \({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.

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

   \({Anti\_CytSec\_MTX}\) is the proportional reduction in the secretion of the proinflammatory cytokines

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

   \({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-{Anti\_CytSec\_MTX})\), (1 + \({Pro\_CytSec\_MTX})\) and \((1-{Anti\_CellInflux\_MTX})\), respectively.

   1. a.

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

      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}/({TCZDrug\_Synovium}+{KD\_TCZ}).\)

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.

1. 7.

   Development and characterization of the Vpop:

   1. a.

      The Vpop was generated using a standardized workflow depicted in Fig. 10.

   2. b.

      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.

      

      The DAS28-CRP score distribution in the Vpop (A) at baseline, (B) post MTX, (C) post Adalimumab and (D) post Tocilizumab therapies. The red dot indicates the mean while the red line indicates standard deviation. The figures show the distribution of DAS28-CRP scores before and after each therapy. The average baseline DAS28-CRP scores were matched to the mean and standard deviation reported in clinical trials (See Supplementary Data 1). Post therapy DAS28-CRP scores indicate that there is a broad distribution of responses seen in the simulations; indicative of phenotypic diversity in the Vpop.

   3. c.

      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.

      

      A Distribution plots for the parameters varied to generate the Vpop (1-64) (B) Distribution plots for the parameters varied to generate the virtual population (65-129). The parameter distribution shows a broad range, indicating that there is sufficient phenotypic diversity in the Vpop. See Table 4 for the key to read these plots and the bounds for the distribution plot corresponding to each parameter.

   4. d.

      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.

      

      Distribution plots for the Baseline Cell Densities in the Vpop for (A) FLS, (B) Endothelial cells, (C)Macrophages, (D) Th1,(E) Th17,(F) CTL, (G) Bcells, (H)Plasma cells and (I)Tregs respectively. The mean/ median values of these cell densities and their ranges are derived from literature (see Supplementary Data 1). The plots show that Vpop has a broad distribution of cell densities at baseline that spans the range determined from the literature. See Table 5 for the ranges for these distributions for cells.

      

      Distribution plots for the Baseline cytokine concentrations in the Vpop for (A)VEGF (B)RANTES (C) TGF-b (D) TNF-a (E) IL-23 (F) IL6 (G)MCP-1 (H) MIP-3a (I)IL-1b (J)GMCSF (K)IFN-g (L)IL-10 (M) IL-12 (N) IL-17 (O)BAFF (P)CAM (Q) AutoAb. The mean/ median values of these cytokine concentrations and their ranges are derived from literature (see Supplementary Data 1). The plots show that the Vpop has a broad distribution of concentrations at baseline that spans the range determined from the literature. See Table 6 for the ranges for these distributions for cytokines.

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.