Abstract
A mechanism-based model was developed to describe the time course of arthritis progression in the rat. Arthritis was induced in male Lewis rats with type II porcine collagen into the base of the tail. Disease progression was monitored by paw swelling, bone mineral density (BMD), body weights, plasma corticosterone (CST) concentrations, and tumor necrosis factor (TNF)-α, interleukin (IL)-1β, IL-6, and glucocorticoid receptor (GR) mRNA expression in paw tissue. Bone mineral density was determined by PIXImus II dual energy X-ray densitometry. Plasma CST was assayed by high-performance liquid chromatography. Cytokine and GR mRNA were determined by quantitative real-time polymerase chain reaction. Disease progression models were constructed from transduction and indirect response models and applied using S-ADAPT software. A delay in the onset of increased paw TNF-α and IL-6 mRNA concentrations was successfully characterized by simple transduction. This rise was closely followed by an up-regulation of GR mRNA and CST concentrations. Paw swelling and body weight responses peaked approximately 21 days after induction, whereas bone mineral density changes were greatest at 23 days after induction. After peak response, the time course in IL-1β, IL-6 mRNA, and paw edema slowly declined toward a disease steady state. Model parameters indicate TNF-α and IL-1β mRNA most significantly induce paw edema, whereas IL-6 mRNA exerted the most influence on BMD. The model for bone mineral density captures rates of turnover of cancellous and cortical bone and the fraction of each in the different regions analyzed. This small systems model integrates and quantitates multiple factors contributing to arthritis in rats.
Models describing the time course of disease progression are becoming increasingly important in drug development because they provide a mathematical structure by which to evaluate drug efficacy, make clinical predictions, and better discern how drugs may be used to target specific markers alone or in combination while avoiding toxicities (Lalonde et al., 2007; Lesko, 2007). Development of these mathematical systems for drug effect analysis is relatively recent. Chan and Holford (2001) reviewed drug treatment effects on disease progression in multiple disease scenarios, placing emphasis on understanding the natural disease progression in the absence of drug. Post and Danhof presented a variety of equations for describing disease progression using pharmacodynamic indirect response models, a mechanistic paradigm for describing turnover of response and drug effect (Dayneka et al., 1993; Post et al., 2005). Other diseases have been described, including bacterial cell growth and irreversible loss after drug therapy, fasting plasma glucose in diabetes, and pain and bone mineral density progression in osteoarthritis (Jusko, 1971; Ravn et al., 1996; Meunier et al., 1999; Frey et al., 2003; Pillai et al., 2004; de Winter et al., 2006). Mechanism-based disease progression models are continually being derived because their insight into how various pathologies contribute to progression and how drugs exert their effects on specific processes is essential to optimizing therapy (Holford and Peace, 1992a,b; Nemeroff, 1994; Hoogendoorn et al., 2005; Holford et al., 2006; Meier et al., 2007; Woo et al., 2008). To date, there are limited mechanism-based models relevant to osteoarthritis and none specific to anti-inflammatory drug effects in rheumatoid arthritis (RA). Furthermore, low-dose corticosteroid therapy has been suggested recently to halt joint erosion in RA (Svensson et al., 2005; Wassenberg et al., 2005; Da Silva et al., 2006). However, before corticosteroid mechanisms can be appreciated in RA, the natural disease baseline must be aptly described by the relevant molecular factors regulating the disease.
This investigation aims at 1) developing a mechanistic model to describe the natural time course of the molecular factors that drive chronic inflammatory arthritis in the rat and that are relevant to corticosteroid effect (this study) and 2) integrating drug binding assumptions into the disease progression model to explain pharmacodynamic effects of corticosteroids (CS) and further resolve disease progression parameters (Earp et al., 2008). The molecular factors assayed to describe disease progression include proinflammatory cytokine mRNA for tumor necrosis factor (TNF)-α, interleukin (IL)-1β, and IL-6 and factors relevant to endogenous corticosteroid dynamics, including glucocorticoid receptor (GR) mRNA and plasma corticosterone (CST). It is well established that TNF-α, IL-1β, and IL-6 increase immune cell trafficking and proliferation and cause up-regulation of GR and plasma CST (Turnbull and Rivier, 1999; Choy and Panayi, 2001; Neeck et al., 2002). In turn, GR and CST are thought to suppress cytokine expression by inhibiting the action of transcription factors nuclear factor-κB and activator protein 1 (Turnbull and Rivier, 1999; Choy and Panayi, 2001; Neeck et al., 2002). These inter-regulation processes between cytokines and endogenous corticosteroid account for limiting inflammation during progression, and they provide common factors for which dexamethasone (DEX) drug effect may be observed and used to explain treatment effects on disease endpoints. Overall disease progression is monitored by edema in the paw, and by femur and lumbar bone mineral density (BMD). The production of edema is directly correlated with immune cell infiltration and proinflammatory cytokine concentrations. Bone production and loss has also been related to the presence of TNF-α, IL-1β, and IL-6 (Rifas, 1999; Gravallese, 2002; Lerner, 2004; Summey and Yosipovitch, 2006; McCormick, 2007). These cytokines not only reduce osteoblast activity but also are thought to increase receptor activator of nuclear factor-κB expression and subsequent activation of osteoclasts.
This study presents one of the first mechanistic models of autoimmune arthritis progression in the rat as described by proinflammatory cytokine mRNA, GR mRNA, CST, paw edema, and BMD. Equations were constructed with the fifth-generation model of corticosteroid dynamics, and they were extended beyond with transduction and indirect response models describing disease initiation and the turnover of primary disease mediators (Ramakrishnan et al., 2002). Model equations were fitted to data, and simulations were performed to identify the contributions of each cytokine to paw edema and BMD. These pathways provide insight into the extent that both corticosteroids and other therapies targeting cytokines may affect reduction in edema and protect BMD during progression of RA.
Materials and Methods
Animals. Male Lewis rats, age 6 to 9 weeks, were purchased from Harlan (Indianapolis, IN), weight matched to approximately 150 g. Animals were housed individually in the University Laboratory Animal Facility and acclimatized for 1 week under constant temperature (22°C), humidity (72%), 12-h light/12-h dark cycle. Rats had free access to rat chow and water. All protocols followed the Principles of Laboratory Animal Care (Institute of Laboratory Animal Resources, 1996) and were approved by the University at Buffalo Institutional Animal Care and Use Committee.
Induction of Collagen-Induced Arthritis in Lewis Rats. The induction of collagen-induced arthritis (CIA) in Lewis rats followed protocols and reagents supplied by Chondrex, Inc. (Redmond, WA). Porcine collagen type II (2 mg/ml) in 0.05 M acetic acid was emulsified with incomplete Freund's adjuvant (IFA; Sigma-Aldrich, St. Louis, MO) using an electric homogenizer (VirTis, Gardiner, NY) equipped with a small blade 10 mm in diameter. Equal amounts of collagen (2 mg/ml) and IFA were mixed in an ice water bath, adding the collagen dropwise to the IFA at the lowest speed setting. The VirTis homogenizer speed was increased to 30,000 rpm for 2.5 min, and then 0 rpm for 2.5 min, and a final mix at 30,000 rpm for 2.5 min. The emulsion was ready when it seemed to be a stiff white substance that congealed instead of dissipating when dropped in water. Ensuring proper time for the solution to cool in the ice bath is critical to prevent collagen degradation (2.5 min was used in between homogenizations). Rats were anesthetized with ketamine/xylazine (75:10 mg/kg) and received 0.2 ml of collagen emulsion by intradermal injection at the base of the tail. Booster injections were administered on the 7th day of the study with 0.1 ml of emulsion to the same injection site.
Experimental Design. On study day 0, 48 rats received injections with collagen emulsion. Animals received booster injections on day 7. Animal body weights and paw edema were monitored on study days 0, 7, 9 to 26, 28, and 30. Three to four rats with arthritis (depending on the availability of rats that had developed arthritis symmetrically in both paws) were killed on study days 9, 12, 15, 17, 19, 21, 23, 25, and 30. Healthy control animals were housed alongside arthritic animals throughout the duration of the study. One healthy rat was sacrificed on study days 9, 15, 19, 23, and 30. Animals were sacrificed by aortal exsanguinations, and blood was collected in syringes containing sufficient EDTA to yield 1.5 mg/ml (4 mM) final concentration (Samtani and Jusko, 2005). Samples were centrifuged at 1800g for 10 min at 4°C. Plasma was collected and stored at –80°C. Skin was removed from hind paws, and the paws were cut off above the ankle and flash-frozen in liquid nitrogen before storage at –80°C.
An additional 20 rats were induced with arthritis for assessment of bone mineral density in the femur, paw, and lumbar regions. Fourteen of the animals that developed the edema in both paws were scanned with the PIXImus II instrument (GE Healthcare, Chalfont St. Giles, UK) to assess changes in bone mineral density. Ten of these rats were scanned for the femur region only on days 8, 10, 13, 15, 19, 22, 27, and 35. Four additional rats were scanned for the femur, lumbar, and paw regions on days 2, 9, 12, 15, 19, 23, and 27. Four healthy controls were also scanned for the femur, lumbar, and paw regions on days 1, 8, 16, 19, 23, 27, and 35.
Paw and body weight data and tissues from preliminary pilot studies in our laboratory and from Earp et al. (2008) were also included in the model development. In total, 162 healthy and arthritic rats were used for all paw mRNA measurements. From this total, 220, 214, and 212 data points were used for TNF-α, IL-1β, and IL-6 mRNA (one to two paws were used per rat). Plasma corticosterone was assayed from 177 rats (177 data points), and paw edema was determined from 314 rats (4340 data points; 13 points per rat). Bone mineral density was assayed from 46 animals, and 1464 data points in total from five different scanned regions were modeled. Only the paw mRNA measures and plasma CST were collected by sparse sampling means. The noninvasive measures, paw edema and bone mineral density, were collected serially.
Measurement of Edema and Body Weight. Edema was indicated by swelling of the rat's hind paws. Two cross-sectional areas were determined with digital calipers (VWR Scientific, Rochester, NY), one area on the forefoot (paw) and the other area at the ankle. Two measurements were made on each section, perpendicular to each other, to define the length and height of the ellipse from which the area was determined. Measurements were made side to side and top to bottom across the paw at the base of the last footpad. Measurements on the ankle were made side to side and front to back at a 45° angle across the ankle. The area contained in the ellipse is as follows: where a is the length of the side-to-side measurement and b is the length of either the front-to-back or top-to-bottom measurements by the calipers. Edema was indicated by the sum of the paw and ankle area measures for each hind foot.
Extraction of Total RNA from Rat Hind Paws. Frozen whole paws (excluding skin) cut off above the ankle were ground with a mortar and pestle under liquid nitrogen. Right and left hind paws were analyzed individually. Total RNA was extracted from the pulverized tissue using column purification with RNeasy mini kits (QIAGEN, Valencia, CA). Manufacturer's instructions were followed with inclusion of an on-column proteinase-K digestion (Ambion, Austin, TX). Tissue samples were homogenized in an ice-water bath for two 15-s durations, with a 20-s period in between to cool. Homogenizations were done with a Kinematica Polytron homogenizer (model PT10-35; Kinematica Inc., Newark, NJ) on speed setting 3.5. After column purification, DNAase digestions were performed with RQ-1 DNAase from Promega (Madison, WI). Sample tubes with enzyme were incubated at 37°C for 30 min. Samples were repurified with the RNeasy columns. The purity and the integrity of the RNA were confirmed by the ratio of optical densities at 260 and 280 nm and by agarose formaldehyde gel electrophoresis. Total RNA concentrations were determined with 2-μl samples in triplicate using a NanoDrop spectrophotometer (model ND-1000; Thermo Fisher Scientific, Wilmington, DE). RNA samples were diluted in nuclease-free water to 25 ng/μl and stored at –80°C before quantitative real-time PCR analysis.
Construction of cRNA Standards for TNF-α, IL-1β, and IL-6. Regions consisting of 708 bp for rat TNF-α (accession no. NM_012675), 807 bp for rat IL-1β (accession no. NM_031512), and 996 bp for rat IL-6 (accession no. NM_012589) were cloned into pCR II-TOPO (Invitrogen, Carlsbad, CA). Primers for TNF-α, IL-1β, and IL-6 were custom synthesized by Integrated DNA Technologies, Inc. (Coralville, IA). Primers for TNF-α cloning were forward 5′-ATGAGCACGGAAAGCATGA-3′ and reverse 5′-TCACAGAGCAATGACTCCAAA-3′. Primers for IL-1β cloning were forward 5′-AAAGAATTCTGAAGCAGCTATGGCAACTG-3′ and reverse 5′-GGCAAGCTTGGGTTCCATGGTGAAGTCAA-3′. Primers for IL-6 cloning were forward 5′-CTGTCTCGAGCCCACCAG-3′ and reverse 5′-TGCAAGAAACCATCTGGCTA-3′. The cRNA standards specific to TNF-α, IL-1β, and IL-6 were synthesized from the cloned constructs by in vitro transcription. The pCR II-TOPO plasmid-clone constructs were linearized using BamH1 followed by in vitro transcription using the MEGASCRIPT T7 polymerase kit (Ambion). Purity of cRNA was checked by the ratio of 260- and 280-nm optical densities and electrophoresis on 5% acrylamide/8 M urea gels. Concentrations of transcribed cRNA clones were determined by absorbance at 260 nm. The GR cRNA standards were cloned previously (Sun et al., 1998).
Quantitative Real-Time PCR. Quantitative real-time PCR was carried out using the Stratagene MX3005P fluorescence-based thermal cycler (Stratagene, La Jolla, CA). TaqMan primers and probes specific to TNF-α, IL-1β, and IL-6 were designed with Primer Express software (Applied Biosystems, Foster City, CA) and custom synthesized by Biosearch Technologies (Novato, CA). Table 1 indicates primer and probe sequences; fluorescent reporter dyes on the 5′ end of the probe; and relevant assay concentrations for probes, primers, and MgCl2. The assay was run with the Brilliant 1-Step Quantitative RT-PCR Core Reagent kit (Stratagene). PCR was conducted with the following cycle parameters: 45°C for 30 min, 95°C for 10 min, and 40 cycles of 95°C for 15 s and 60°C for 1 min. Standard curves were generated for each cytokine mRNA assay with cloned cRNA standards, which were assayed in duplicate for each experiment and samples were assayed in triplicate. A single additional control tube was run for each sample excluding the reverse transcriptase to test for genomic DNA contamination. No significant amplification was observed for all control tubes. Inter- and intraday assay variation was less than 15% coefficient of variation for all quantitative real-time PCR assays.
Corticosterone Extraction and High-Performance Liquid Chromatography Analysis. Rat plasma was thawed, and 0.1- to 1.0-ml aliquots were extracted with methylene chloride (MeCl2) in 7-ml glass Pyrex tubes (Corning Inc., Corning, NY). Tubes were shaken for 45 min before washing the MeCl2 phase with 0.5 ml of 0.1 N sodium hydroxide. The NaOH phase was removed after centrifugation, and the MeCl2 phase was washed twice with 0.5 ml of water, discarding the water after each wash. MeCl2 was evaporated off with purified air, leaving a residue that was reconstituted in 110 μl of mobile phase. Chromatography conditions involved a mobile phase of 600 ml of MeCl2, 350 ml of heptane, 10 ml of glacial acetic acid, and 54 ml of ethanol; a Zorbax normal phase column (Agilent Technologies, Santa Clara, CA); Waters model 1515 isocratic pump; and Waters model 2487 dual-wavelength absorbance detector (Waters, Milford, MA) (Haughey and Jusko, 1988). The lower limit of quantification was 5 ng/ml with an intraday coefficient of variation of less than 10%.
Measurement of Bone Mineral Density. BMD measurements were performed using the PIXImus II dual energy X-ray absorptiometer (GE Healthcare). Rats were anesthetized with 4% isoflurane and maintained at 1.5% isoflurane for the duration of the scan. Animals awakened within 30 s to 1 min after removal from isoflurane. Before the investigation started an instrument calibration was run. Quality controls were examined daily using the supplied phantom (BMD = 0.0664 g/cm2; percentage of fat = 10.4%). Anesthetized rats were positioned on the PIXImus trays, and their tail or hind paw was taped to secure the rat in three positions for determination of BMD in the femur, lumbar, and paw regions of interest.
The 90° position was used for the femur scans as proposed by Binkley et al. (2003) and applied by Soon et al. (2006). The position for the lumbar region of interest was also adapted from these two references. The paw scans were done with the rat lying on its side, back away from the machine, and the upper leg taped back so the bottom leg could rest freely with the paw on its side. Scans of the whole paw and lower leg from the side were obtained. Paw scans were included to determine whether the localized edema in the paw affected BMD more than in the femur or lumbar regions. BMD was quantified by the Lunar PIXImus II software version 2.00, with a threshold value of 1345. Four regions of interest in the femur were analyzed, and they are shown in Fig. 1, with the corresponding pixel dimensions of each region of interest.
Mechanism-Based Disease Progression Model.Figure 2 shows a general schematic for the entire arthritis progression model. All compartments in this model are described by differential equations that maintain homeostasis when unperturbed. Disease progression is initiated by the transit compartments preceding each cytokine mRNA.
Corticosterone and GR dynamics. The fifth-generation model for corticosteroid dynamics proposed by Ramakrishnan et al. (2002) and modified by Hazra et al. (2007a,b, 2008) was adapted and applied to the effects of endogenous corticosterone via the glucocorticoid receptor on cytokine mRNA production, GR mRNA production, paw edema, and bone mineral density. The equation describing GR mRNA turnover is as follows: where kin_GRm is the first-order production rate of GR mRNA. Drug bound to receptor is indicated by DR. After the DR complex translocates to the nucleus it is defined as DRN. The first-order translocation constant (kT = 58.2 h–1) between DR and DRN was derived previously (Ramakrishnan et al., 2002; Hazra et al., 2007a,b, 2008) and is fixed for all CS in this model. IC50_GRm/DRN is the concentration of DRN required to inhibit GR mRNA production by 50%, and γDR is the Hill coefficient for inhibition by DRN. The rate constant kout_GRm describes the first-order loss for GR mRNA and is defined through the steady-state baseline condition: The variable Tcyt is representative of an overall cytokine protein expression and accounts for the slow onset of GR mRNA relative to the rapid rise in cytokine response. The differential equation for Tcyt is denoted by the following: where Acyt is the contribution of proinflammatory cytokine mRNA on the production of Tcyt: The α1 and α2 are the intrinsic activities of each cytokine mRNA for producing cytokine protein, Tcyt. Correlation with IL-6 cytokine mRNA in previous model generations indicated that IL-6 contribution to GR mRNA was negligible (i.e., α3 for IL-6 < 10–8). This was fixed to zero for simplicity in this model.
Depletion of the free GR due to translocation to the nucleus, degradation, limited synthesis rate, and redistribution play major roles in governing the shape of time course of effects from corticosterone and synthetic glucocorticoids. The equations describing concentrations of free GR are as follows: where ksyn_GR is the first-order synthesis rate from GR mRNA and kdgr_GR is the first-order loss rate constant. DRN_C is the amount of CST-GR complex in the nucleus, and DRN equals DRN_C for the disease progression model without drug effect. The parameter kRE_C describes the rate of which corticosterone bound to receptor returns from the nucleus, and RF is the fraction that returns intact and active. The values for kdgr_GR and RF were fixed from the literature (Ramakrishnan et al., 2002; Hazra et al., 2007a,b, 2008). The second-order binding rate constant for CST was defined as follows: where 346.47 ng is the molecular mass of CST and 18.64 nM is the literature reported KD value for CST (Buchwald, 2008). The parameter ffGC · kon_C accounts for the fraction of drug in plasma that is capable of binding in tissue. The model assumes that CST equilibrium between tissue and vasculature is instantaneous. Bound CST-receptor complex in the cytosol is defined as follows: where kT is the first-order rate constant for translocation of bound CST-GR complex to the nucleus as reported previously (Hazra et al., 2007a,b, 2008). Concentrations of bound CST-GR complex in the nucleus are described by the following: This notation is important for pharmacodynamic models when other corticosteroids are present. In this case, when no drug is present: DRN = DRN_C. The parameters GR0, DR0, and DRN0 are the amounts of the glucocorticoid receptor, corticosteroid bound to receptor, and corticosteroid bound to receptor in the nucleus at time 0 under baseline conditions at steady state and when kdis values are set to 1. These baseline values at time 0 are defined as follows:
Plasma CST concentrations are described by the following: where kin_CST is the first-order rate constant for the production of CST, and BCYT is the influence of proinflammatory cytokines on CST up-regulation defined as where β1 and β2 are the intrinsic activities of each individual cytokine on CST up-regulation. Correlation with IL-6 cytokine mRNA in previous model generations indicated that IL-6 contribution to plasma corticosterone was also exceedingly small (i.e., β3 for IL-6 < 10–8) given the linear model structure and β3 for IL-6 was fixed to zero for simplicity in this model.
Proinflammatory cytokine dynamics. For animals that are induced with the disease, the production rate (kt · R0) for cytokine mRNA is increased to a new disease steady-state value (kt · R0_mRNA · kdis_mRNA) at time of induction, time = 0 h. For the natural healthy state, the proinflammatory cytokines are never turned on and kdis_mRNA maintains a value of 1.0. The equations and initial conditions describing the first event compartment for cytokine mRNA regulation for time ≥0 h are as follows: where T1 indicates the first transit compartment relevant for the series of events that takes place before when disease onset is observed. The values of kt indicate the turnover rate constant for each compartment in this series. Baseline values for the cytokine mRNAs are indicated for when time <0 h and kdis_mRNA is equal to 1. The final transit compartments in each series are defined by the following: where subscripts denote the relevant mRNA. Transduction models with different numbers of compartments were tested for each cytokine mRNA independently of each other or other model components to capture the number of compartments required to account for both a large delay and rapid rise to a new disease steady state using data from animals without DEX treatment. The meaning of each compartment is arbitrary other than being a component of a series of events that take place during the onset of arthritis in the rats. This interpretation is more true to the underlying physiology than a delay or onset time that is also more numerically difficult to compute. Onset times additionally do not provide the smooth rise observed early in the onset of cytokine mRNA expression.
Cytokine mRNA is governed by many different factors, and thus equations describing the production and loss of each cytokine are different. The turnover of TNF-α mRNA is described by the following: where kin_TNF describes the first-order production rate dependent upon transit compartment T29TNFα. DRN is endogenous corticosterone bound with receptor in the nucleus, IC50_TNFα/DRN is the concentration of DRN required to inhibit TNF-α production 50%, and γ1 is the Hill coefficient for the inhibition by DRN. The first subscript for each IC50 parameter after the underscore denotes the factor being altered. The subscript for each IC50 parameter following the back-slash indicates the factor that is causing the inhibition. The term IL6mRNA describes the amount of IL-6 mRNA, and IC50_TNFα/IL6 denotes the amount of IL-6 mRNA necessary to inhibit TNF-α production 50% of maximal capacity (0.3). The rate constant kout_TNFα describes the first-order loss for TNF-α mRNA and is defined through the steady-state baseline condition as follows: where DRN0 and IL60 describe the baseline values of DRN and IL-6 mRNA.
The turnover of IL-1β mRNA is described by the following: where kin_IL1β describes the first-order production rate dependent upon transit compartment T19IL1β. The IC50_IL1β is the concentration of DRN required to inhibit IL-1β mRNA production 50%. The term Rem describes an increase in other factors that control the decline in inflammation and is defined by Rem = T27IL1β – R0_IL1β. The parameter Rem was important to account for an apparent remission in IL-6 and IL-1β mRNA that CST alone could not explain. Changes in IL-6 and IL-1β mRNA resulting from high CST concentrations were not as evident as those from DEX. The apparent lack of sensitivity of these factors to CST was important in determining how much effect CST had on the biomarkers after DEX was administered and CST concentrations plummeted. If CST entirely explained the degree of remission for both cytokines, then in Earp et al. (2008) after DEX was administered predicted IL-6 and IL-1β mRNA concentrations would spike well above the observed baseline when DEX concentrations fell and CST concentrations remained low. The parameter IC50_IL1β/Rem denotes the amount of Rem necessary to inhibit IL-1β mRNA production to 50% of maximal capacity (0.7). The rate constant kout_IL1β describes the first-order loss rate-constant for IL-1β mRNA and is defined through the steady-state baseline condition as follows: The turnover of IL-6 mRNA is defined by the following: where kin_TNFα describes the first-order production rate dependent upon transit compartment T24IL6. The parameter IC50_IL6/DRN is the concentration of DRN required to inhibit IL-6 mRNA production by 50%, and γ3 is the Hill coefficient for inhibition by DRN. The rate constant kout_IL6 describes the first-order loss for IL-6 mRNA, and it is defined through the steady-state baseline condition as follows:
Dynamics of paw edema. Paw edema was modeled as production from natural growth (kgrowth), production from cytokines (kin_Paw), and loss of the edema (kin_Paw/R0_Paw) that was produced by cytokines (i.e., when edema is not present only kgrowth affects paw size) as follows:
Pcyt describes the total contribution of cytokine mRNA to the overall edema. The values of π for each cytokine reflect intrinsic activities of each cytokine on the production of paw edema.
Dynamics of bone mineral density. Bone mineral density was modeled from five different regions of interest (four femur, one lumbar). The major assumptions governing this model are 1) there are two types of bone, cancellous (Canc) and cortical (Cort), that contribute to BMD in each scan and 2) each region analyzed is composed of different fractions of cancellous and cortical bone. The parameters from the equations for the two types of bone turnover are the same for all regions; the only difference is the fraction of cancellous bone that is modeled for each specific region. Equations for cancellous and cortical bone mineral density are as follows: and they are based on the logistic function with modifications. The parameters kgrowth and BMDmax indicate the first-order growth constant and maximum value of bone density for each type. Two hypothetical compartments were created to account for effects of cytokines and corticosteroids on BMD. Proinflammatory cytokines are thought to exert their effects on bone loss primarily through the up-regulation of receptor activator of nuclear factor-κB and activation of osteoclasts. Hypothetical osteoclast (OC) activity was described by the following compartment with a baseline of 1.0 for simplicity. where Rcyt is defined as follows: and each ρ value indicates the individual contribution of that cytokine to the production of OC, which stimulates bone loss in eqs. 27 and 28.
A hypothetical compartment for osteoblasts (OB) was assumed that is thought to be reduced in number and activity by corticosteroids. The effect of DRN on BMD was delayed by nine transit compartments before inhibiting osteoblast production. The baseline was fixed to 1.0 for simplicity. Model equations are as follows: where kOB describes the turnover rate of the hypothetical osteoblast compartment. The value of IC50_GC is the concentration of DRN required for 50% inhibition of “osteoblast” production.
Although the above-mentioned equations were used to model the turnover of cortical and cancellous bone, BMD in each region was modeled as a function of cancellous and cortical bone densities and the fraction of bone in each region that is cancellous (FTF, FDF, FMF, FEF, FLR,) and the fraction of bone that is cortical (1 – Fi). Equations for each specific region are as follows: Subscripts are indicative of each specific region (TF, total femur; DF, diaphyseal femur; MF, metaphyseal femur; EF, epiphyseal femur, and LR, lumbar region) or bone type (cancellous or cortical).
Model Fitting and Nonlinear Regression Analysis. Owing to the model complexity and amount of paw edema and BMD data compared with cytokine data, the model was constructed in two phases: 1) fitting of molecular biomarkers and 2) fitting disease endpoints. In addition, all dexamethasone drug effect data presented in Earp et al. (2008) were incorporated and fitted in the disease progression analysis simultaneously. The first segment of the model to be fitted were disease and drug effects on cytokine mRNA, GR mRNA, and plasma corticosterone. The relevant parameters were obtained and then fixed for the model fitting of paw edema and BMD data. Arthritis disease progression-related data, parameters, and model-fitted curves are presented here.
Model fittings were performed using the maximum likelihood algorithm in S-ADAPT (Biomedical Simulations Resource, University of Southern California, San Diego, CA) with the OPTIMIZE module. The pooled fit command was used to permit a NONMEM type data set to be used with the lsoda algorithm for solving differential equations. The Nelder-Mead simplex search algorithm was implemented for the model relevant to the molecular biomarkers and paw edema because this was an exploratory analysis and initial estimates are difficult given the lack of literature data for these biomarkers and the complexity of the model. The Broyden-Fletcher-Goldfarb-Shanno variable metric method was used for the BMD model because convergence is achieved faster than the Nelder-Mead simplex and initial estimates were more robustly obtained from features of the BMD data and from literature values (Bilezikian et al., 2002). Both algorithms converged successfully. Simulations and model code from S-ADAPT were cross-checked using Berkeley Madonna Software (University of California, Berkeley, CA), and all simulation data for figures were done with Berkeley Madonna implementing the Rosenbrock (stiff) method for solving differential equations. All data were pooled, and initial attempts at a population analysis using the complete model for all measured biomarker responses were not successful due to sparse and variable mRNA data and computational limitations owing to the size and complexity of the model.
Results
Disease Progression of TNF-α, IL-1β, and IL-6 mRNA. The disease progression for proinflammatory cytokine mRNA for TNF-α, IL-1β, and IL-6 is shown over 960 h (34 days) after induction of CIA in Lewis rats in Fig. 3A. Marked delays of 240, 288, and 360 h were observed before rapid rises in IL-6, TNF-α, and IL-1β mRNA profiles. Expression of TNF-α and IL-1β mRNA expression rose to plateaus of 110 and 97 ng mRNA/mg total RNA, and remained elevated for the remainder of the monitored time course. In contrast, IL-6 mRNA rose to an abrupt peak at 6225 ng mRNA/mg total RNA on day 19 (456 h) before declining to an apparent steady-state disease amount of approximately 1300 ng mRNA/mg total RNA. Interestingly the expression of IL-6 mRNA was induced to a much larger extent than either TNF-α or IL-1β mRNA, peaking at 6225 ng mRNA/mg total RNA at 456 h after induction. The interanimal percentage coefficients of variation were between 35 and 60% for each time point in the arthritic groups, with five to nine samples per time point. Figure 3B shows the baseline mRNA expression in healthy rats. Baseline values for cytokine mRNAs were determined as the average of the observed values for all time points from healthy rats, and they are presented in Table 2. The number of available healthy samples at each time point was fewer at n = 2 to 4. No time dependence in expression of cytokine mRNA from healthy rats was observed.
Model fittings in Fig. 3 seem to fit the data well. The number of transit compartments necessary to account for the disease progression varied between cytokine type. TNF-α mRNA required 29 compartments to account for the delay, with a faster transit time to account for the abrupt rise. The rise in IL-1β mRNA expression was slower and delayed further than TNF-α or IL-6; 19 compartments were necessary to describe the onset of IL-1β mRNA, with a longer transit time. The onset of IL-6 mRNA was the earliest, with the fastest rise of all the measured cytokine markers, and 24 transit compartments were needed. Parameter estimates for cytokine disease progression are shown in Table 2. Baseline values for each cytokine were fixed to the average of all measured responses from each time point in healthy rats. The IC50 values indicated the sensitivity of the cytokine to CST-GR complex in the nucleus. Interleukin-6 mRNA seemed to be the most sensitive to DRN (IC50_IL6 = 4.50 nM). The IC50 value for IL-1β was 7-fold higher at 32.5 nM and even higher for TNF-α at 550 nM. The kdis values indicate the maximum mRNA expression in the disease state without the presence of endogenous CST and DRN inhibiting the production of these cytokines.
GR mRNA and Plasma CST Concentrations.Figure 4 shows the disease progression time course for GR mRNA (Fig. 4A) in paw tissue and plasma corticosterone concentrations (Fig. 4B). In both cases, a delay of 240 h was observed before the onset of response to CIA. The GR mRNA rose slowly, and it did not reach a plateau during the monitored time course of disease progression. The highest observed response was the second-to-last measured time point, with a mean value of 620 ng mRNA/mg total RNA. CST rose rapidly at the same time with a similar trend to TNF-α mRNA. Mean CST concentrations did not exceed 350 ng/ml. Model parameters for the fifth-generation model of corticosteroid dynamics and for the turnover of GR mRNA and CST are presented in Table 3. Baseline values were fixed to the average of the data from healthy animals. Intrinsic activity parameters for effect of cytokine mRNA on GR mRNA or plasma CST showed that TNF-α had the greatest effect on these two responses during disease progression. The value of α3 and β3 for IL-6 was set to zero because there was a decline in response of IL-6 mRNA as the disease progressed, and this decline was not observed for GR mRNA or CST concentrations. This is probably why estimates of α3 and β3 tended toward zero (estimates resulting in values less than 10–8) during early model fittings. Although it is known that high CST concentrations suppress GR mRNA, this was not readily observed in the data because a continual rise in GR occurred, despite a continual rise in plasma CST. The fact that GR concentrations increase while cytokine concentrations reach a plateau suggests that GR mRNA increases, in part, because of genomic regulation and not simply because of immune cell infiltration into paw tissue. In addition, the IC50 value for CST effects on GR mRNA was estimated at 550 nM, similar to the value for TNF-α mRNA inhibition. Model fittings are shown as lines in Fig. 4 and seemed to capture the disease trend well.
Paw Edema Progression. The disease progression of CIA as indicated by paw edema is shown in Fig. 5. Similar to cytokine mRNA expression, there is a delay in onset of 288 h, a peak of 140 mm2, and slight remission to a response of 120 mm2 during CIA progression. The increase in mean observed paw size from natural animal growth was modest, increasing from 77 to 86 mm2. The model was described by a first-order production and first-order loss of edema dependent upon the presence of proinflammatory cytokines and the extent of edema, and it successfully described the data. Model parameters are shown in Table 4. Values of π1, π2, and π3 indicate the extent to which TNF-α, IL-1β, and IL-6 mRNA contribute to the edema. Larger π values for TNF-α and IL-1β mean less mRNA is required to contribute to the edema than for IL-6 where more expression is required for the same response, similar to a measure of potency. Being dependent on cytokine mRNA production, the remission observed in paw swelling is then explained by the remission observed in IL-6 mRNA expression. The CV% values for paw parameters were low but dependent upon the validity of the fixed cytokine mRNA model.
Simulations of the paw edema model were done with and without inhibition of different cytokine effects on edema to indicate how each cytokine mRNA contributes to the overall edema, relevant to newer RA therapies that target blocking of cytokine effects. Figure 6 shows these simulations. Figure 6, A to C, shows paw edema response after continual inhibition of TNF-α, IL-1β, or IL-6 mRNA by 0, 50, and 100%. Figure 6, C and D, depicts the paw edema progression when mRNAs of two cytokines are inhibited jointly by 0, 50, and 100%. The 100% joint inhibition of two cytokines yields a profile of paw edema that is indicative of the contribution the uninhibited cytokine has on edema production. These 100% dual-inhibition profiles in Fig. 6, C and D, seem to have the same trend in progression as the uninhibited cytokine due to the linear relationship between mRNA and production of edema.
Bone Mineral Density Progression.Figure 7 shows femoral and lumbar BMD in healthy controls and during CIA in Lewis rats. Scans of ankle and paw BMD yielded no change and high variation, most likely due to multiple overlapping bones when scanned from the side. However, changes in BMD were quite apparent in the femur. Skeletal growth was observed throughout the time course of healthy rats. In arthritic animals, reduced BMD occurred beginning at 420 h after induction, later than the onset of cytokine mRNA. Changes in BMD were less apparent in the lumbar region between the two groups.
In general the model fitted the data well, considering that there were two fairly restrictive assumptions 1) only cancellous and cortical bone contributes to BMD, and 2) differences in each region of interest are solely due to different fractions of cancellous and cortical bone. Parameter estimates for bone turnover are presented in Table 5. Initial values at time 0 and maximal bone densities for each bone type were as expected, being lower for cancellous bone and higher for cortical bone. The parameter estimates for fractions of cancellous bone are reported in Table 6. These estimates were surprisingly close to literature-reported values (Bilezikian et al., 2002) for sections of long bone, with 20 and 30% of bone being cancellous, and values for the distal end of the femur in the knee joint (epiphyseal femur) and the lumbar region were higher at 63 and 70%.
Discussion
Development of a mechanistic model for rheumatoid arthritis is complex owing to the many molecular factors that affect disease progression. The model was derived based on biomarkers that are well known as major mediators of inflammation and bone erosion. It is well established that there is an endocrine link between proinflammatory cytokine production and glucocorticoid receptor and corticosterone regulation. This partly explains disease remission, with increased expression of GR mRNA and CST after CIA onset, and it permits a mechanistic framework to examine the effects of other corticosteroids.
The structure of the disease progression model integrates endogenous glucocorticoid effects with cytokine mRNA expression in an effort to model the effects of DEX in an animal model of arthritis (Earp et al., 2008). DEX is a synthetic glucocorticoid that perturbs many inflammatory pathways and all measured biomarkers acting through the same mechanisms as CST. The turnover of these markers within the disease progression and their effects on disease endpoints can be obtained from the PD data after DEX administration. The disease progression data alone cannot resolve all parameters. Rationalization of these remaining equations and parameters requires information from DEX PD, and it is presented here to a limited extent and in more detail in Earp et al. (2008).
The model describing disease progression in paw cytokine mRNA starts as a simple transduction model with a rate constant to account for the delay, baseline steady-state initial condition for the healthy response, and the kdis rate constant to increase the cytokine mRNA to a new disease steady state at time of induction. This is a new application of transduction compartments to describe four aspects of disease progression: time of induction, delay in disease onset, rate of progression, and new disease steady-state response. Equations for individual cytokine mRNAs were complicated by inhibition of mRNA production by CST-GR complex in the nucleus, possible disease remission, and cross-talk between IL-6 and TNF-α mRNA. Since DEX and CST are assumed to act similarly, many of these modeling complexities were resolved from the DEX PD data (Earp et al., 2008).
Initially the decline in IL-6 paw mRNA was explained by CST binding to GR and DRN inhibiting IL-6 mRNA production as DRN rose above the IC50 value for IL-6 (4.5 nM) mRNA inhibition during disease progression. However, early simulations indicated that CST alone could not account for the entire remission of IL-6. If DRN rose rapidly as in the case of CS administration CST would be suppressed rapidly and for a prolonged time (Earp et al., 2008), such that when DRN returned to normal and CST remained suppressed, the IL-6 mRNA production would overshoot the measured response.
Abrupt DEX PD changes in cytokine mRNA indicated that rate constants relevant to the loss of mRNA were much faster than the transduction rate constants describing disease progression. This observation forced the production and loss constants for cytokine mRNA to be modeled separately from the transduction rate constants. Literature indicates rapid production and loss of cytokines within hours, whereas disease effects may be drawn out over days. Using transduction equations to describe disease progression provides a conceptual link between slower rate-limiting disease steps and measured responses with faster turnover.
The localized nature of the disease in paw/joint tissue made measurements of TNF-α, IL-1β, or IL-6 protein in paw tissue impractical. We tried measuring these cytokine proteins in homogenized paw tissue diluted with phosphate-buffered saline and plasma. Background/matrix effects were a problem, and the TNF-α, IL-1β, and IL-6 enzyme-linked immunosorbent assay kits (Promega) were not valid for tissue analysis. The cytokine enzyme-linked immunosorbent assay kits for protein measurements were not sensitive enough for plasma cytokine concentrations in arthritic rats. Given these experimental limitations, the model uses the available methods to relate physiologically relevant events. Owing to the rapid turnover of both cytokine protein and mRNA in tissue and to these proteins being constitutively expressed in response to an antigenic stimulus, time courses of cytokine mRNA may be reflective of synthesized cytokine. Thus, correlations between cytokine mRNA and effects on disease progression may be useful if the turnover rates of disease endpoints are slower than those of cytokine protein.
The DEX PD data helped resolve Hill coefficients for DRN effect from different doses. For TNF-α and IL-6 mRNA, these values were fixed to 2.0 because this agreed with both the initial model estimations and the physiological homodimerization of CST-GR complex in the nucleus.
Paw edema was monitored in every animal that developed the disease. There were 4943 paw edema data points used in the disease progression, healthy animal profiles, and pharmacodynamic analysis compared with 439 points for mRNA of each cytokine. Due to variation in paw edema progression among animals, describing paw swelling during simultaneous fitting of all data initially biased the model fitting because it placed most emphasis on the how the paw edema was modeled and little on the cytokine mRNA driving the edema. Without capturing the mRNA responses first, the development of the structural model describing cytokine mRNA and paw edema was difficult and unclear. The numerous feedback loops by DRN on most inflammatory factors and these factors on GR and CST meant that if one variable was altered, model descriptions for all other biomarkers were also altered. Thus, the model describing cytokine mRNA, GR mRNA, and plasma CST was developed first and fixed for evaluating the structural model for paw edema dependent on the cytokine mRNA expression. This permitted shorter model fitting times and better model revisions for the complex feedback model between cytokine and CST-GR regulation.
The final model for paw edema was relatively simple and driven by contributions from all three cytokines. Paw edema increased early as affected by TNF-α and IL-6, and the decline in IL-6 mRNA after day 21 helped explain the similar decline to a steady state in paw edema. These individual contributions of cytokines on edema were further resolved by administration of DEX, which inhibited each cytokine mRNA profile at different rates and extents while correlating these mRNA amounts to paw edema.
There is increasing interest for treating rheumatoid arthritis with biologics aimed at inhibiting the effects of TNF-α and IL-1β. Model simulations continually inhibiting the effects of each cytokine by 50 or 100% were performed to predict extent of edema reduction. Because TNF-α, IL-1β, IL-6, and other proinflammatory cytokines have overlapping functions in producing edema, completely suppressing the inflammatory response by targeting only one or two of these pathways results in incomplete reduction of edema. The myriad pathways contributing to inflammation emphasizes the role of corticosteroids in treatment of RA for their inherent ability to suppress many different inflammatory processes simultaneously. Understanding the extent to which dexamethasone, prednisolone, or other CS reduce each inflammatory signaling pathway in light of how other therapies affect these common disease factors may contribute to determining optimal treatment regimens for agents given alone or in combination.
Bone mineral density was modeled in each region by the sum of the cancellous bone plus the cortical bone. Because these amounts were not known, the fraction of cancellous bone was modeled in each region. Physiological understanding that bone turnover is more rapid in cancellous bone is in agreement with model estimates for initial and maximum BMD in both bone types. However, the turnover rate constant for cancellous bone was slightly lower than that for cortical bone. This is caused in part by the parameter ftc that reduces the effect of the OC on cortical bone turnover.
The model does not consider proximity to inflamed joints as a factor, and as such the lumbar region is treated as if equal to the knee joint in terms of inflammatory effects. Although disease effects may not be readily apparent in the lumbar region of the rat, differences are observed in human RA in the spine. The time course of response may need to be carried out longer to observe these inflammatory effects on BMD in the rat lumbar region.
High CV% values for bone parameters were rather surprising because these parameters seemed to stabilize quickly and shifted very little during model convergence. It is possible that the high CV% came from interanimal variation in onset of disease progression and rate of bone loss. Precision measurements for total BMD have been reported previously (Binkley et al., 2003; Soon et al., 2006). In all cases, the CV% value was less than 10%. It is possible this CV% value also comes from not having measured cancellous or cortical bone directly. However, the value of the parameter estimates seemed reasonable.
Trends in both disease progression and pharmacodynamic responses to dexamethasone (Earp et al., 2008) were necessary to discern model behavior of proinflammatory cytokine mRNA and their effect on driving disease progression in paw edema as well as endogenous CST and GR mRNA that in turn regulate the production of TNF-α, IL-1β, and IL-6 mRNA. Other components such as immune cell presence in tissue, amounts of transcription factors in the nucleus, and additional inflammatory regulators (e.g., prostaglandins, nitric oxide, T-helper 2 cytokines) are valuable components to understanding inflammation. Despite the absence of these other factors in the model, edema and BMD seem to be well described by the major cytokines that drive the production of inflammation. This is the first detailed mechanistic model for CIA progression in the rat and model qualification may come with future studies. Defining and developing translational assumptions between disease progression in the rat and clinical disease observations may not be possible for all biomarkers. However, it is possible that for BMD, the extent of disease or drug effect will correlate well, even if the rates of turnover are different because these are chronic conditions with chronic therapies, and disease progression and drug responses tend toward steady-state values over time. Owing to the endogenous link between CST and cytokine mRNA, treatments with DEX and other corticosteroids in the clinic along with predetermined knowledge of drug receptor binding may aid in relevant clinical predictions. Future experiments aim to validate the model by predicting the response to methylprednisolone dosing to arthritic rats. Experiments with other therapeutics may prove informative if the drug effects are limited individual pathways in the model. Anti-cytokine antibodies are a good example of therapeutic agents that may block specific signaling components to identify the relative contributions of certain cytokines on disease endpoints.
Footnotes
-
This work was supported by National Institutes of Health Grant GM24211 and by a fellowship from Amgen Inc. (to J.C.E.).
-
Article, publication date, and citation information can be found at http://jpet.aspetjournals.org.
-
doi:10.1124/jpet.108.137372.
-
ABBREVIATIONS: RA, rheumatoid arthritis; CS, corticosteroid; TNF, tumor necrosis factor; IL, interleukin; GR, glucocorticoid receptor; CST, corticosterone; DEX, dexamethasone; BMD, bone mineral density; CIA, collagen-induced arthritis; IFA, incomplete Freund's adjuvant; PCR, polymerase chain reaction; bp, base pair(s); FAM, 5-carboxyfluorescein; DR, drug bound to receptor; DRN, drug bound to receptor in the nucleus; Canc, cancellous; Cort, cortical; OC, osteoclast; OB, osteoblast(s); TF, total femur; DF, diaphyseal femur; MF, metaphyseal femur; EF, epiphyseal femur; LR, lumbar region; CV%, percentage coefficient of variation; PD, pharmacodynamic.
- Received February 1, 2008.
- Accepted April 10, 2008.
- The American Society for Pharmacology and Experimental Therapeutics