Design and study population
This was an observational, longitudinal study of patients with dcSSc registered in the EUSTAR database. The EUSTAR network has been described elsewhere [14, 15]. In brief, EUSTAR is a growing database of patients with scleroderma treated at centers worldwide; all patients undergo annual scheduled clinic visits, providing observational, longitudinal data.
For this analysis, patient data were extracted from January 1995 to February 2019. Ethics Committee approval was obtained from all centers, and informed consent was provided when required by the ethical regulations at the specific centers. Patients were ≥ 18 years of age, and to have SSc as classified by the American College of Rheumatology (ACR) (1980) or ACR/European League Against Rheumatism criteria (2013) [16, 17]. Patients with dcSSc were identified from this cohort based on the LeRoy criteria , and those with the available data for the first non-Raynaud’s manifestation were extracted. If data on the LeRoy criteria were not available, then the extent of skin involvement (e.g., skin fibrosis at any time with mRSS ≥ 1 of upper arms, thorax, abdomen or thighs) was used as a surrogate.
Population and outcomes
Analyses were performed using patients with dcSSc who met the above criteria with or without HAQ-DI score assessments. Patients with HAQ-DI score were further divided into either those who had ≥ 1 HAQ-DI and those who had ≥ 2 HAQ-DI score assessments. Clinical data used in analyses are shown in Table 1, Fig. 1 and Fig. S1.
Analyses of the impact of organ involvement focused on the major advanced complications seen in dcSSc patients. Advanced gastrointestinal (GI) events were defined as malabsorption or ³ 10% weight loss from baseline. Echocardiographic measurement of systolic pulmonary arterial pressure > 45 mmHg was used as a proxy for pulmonary hypertension (PH) as data on right heart catheterization were not collected routinely in the database. Interstitial lung disease (ILD) was confirmed by high-resolution computed tomography or chest X-ray, and significant lung involvement was defined as ILD on imaging and forced vital capacity < 75% predicted. Cardiac involvement was defined as left ventricular ejection fraction (LVEF) £ 45% measured on echocardiography, and renal involvement was defined as the presence of renal crisis (abrupt onset of severe hypertension accompanied by rapidly progressive renal failure, hypertensive encephalopathy, congestive heart failure, and/or microangiopathic hemolytic anemia).
Patients could be enrolled at any time in the disease course if they were classified as having SSc (defined as baseline). Clinical history, demographic characteristics, use of immunomodulatory treatment, and death were documented for patients with or without HAQ-DI score assessments; this information was based on medical records in the EUSTAR database. It was intended that patients should be documented once a year. Death was documented with the date or year (if the date was not exactly known) and reason of death.
Regression analyses to identify predictors of death and HAQ-DI score progression
Patients with ≥ 2 HAQ-DI score assessments were selected to model prediction of HAQ-DI score at 1 year (i.e., progression). Patients with ≥ 2 specific organ assessments were selected to predict major advanced organ involvement. Using these populations, the relationship between baseline characteristics and death rates and HAQ-DI score progression at 1 year were examined both at the univariate level and in multivariable models. Death rates and HAQ-DI scores were analyzed by Cox regression and linear regression analyses respectively, in relation to each baseline covariate: gender, age at onset of onset of Raynaud’s phenomenon, RNA polymerase III positive, immunomodulator treatment, corticosteroid use > 10 mg/day, HAQ-DI score, mRSS, and number of major advanced organs involved (0, 1, or ≥ 2 based on lung, cardiac, GI, PH, or renal involvement). All analyses were performed using IBM® SPSS statistics 24.0 and R software (version 3.4.4).
Transition model for dcSSc
A Markov model was developed to determine the natural progression of dcSSc over a patient’s lifetime horizon. Within the above-mentioned dataset, three risk equations were developed to model the evolution of the following endpoints: HAQ-DI score; major advanced organ involvement; death.
The HAQ-DI score is a continuous variable, and was transformed into a categoric variable in this model to represent five different health states ([0─0.5], [0.5─1.0], [1.0─1.5], [1.5─2.0], and [2.0─3.0]); the health states were based on expert opinion and were initially defined as uniformly splitting the 0─3 HAQ-DI score range into six categories (i.e., this is the maximum amount of categories that enable the multi-state modelling (MSM) package to converge when fitting a transition matrix between HAQ-DI states); however, the two last categories did not contain enough patients and were grouped into a single state, leaving five states in total. The HAQ-DI transition matrix was calibrated using an MSM that provides constant transition probabilities over time from any HAQ-DI state to any other one. For calibration of the HAQ-DI transition matrix, transition intensities were developed using longitudinal, patient-level data for HAQ-DI states as a function of gender, age at baseline, and lung status at baseline. A 1-year cycle was chosen to transform intensities into yearly transition probabilities of moving between the five HAQ-DI states (Table S1).
For calibration of the organ equations, longitudinal patient-level data were used to calculate survival in the no-organ-involved state in relation to specific patient characteristics at baseline (age, sex, HAQ-DI states, and major advanced organ involvement such as lung, PH, cardiac, renal, or GI) (Table S2). This equation was calibrated using the survival package (R software).
For calibration of the third equation, death rates were calculated using a standardized mortality ratio (SMR), defined as the ratio of observed deaths in the dcSSc EUSTAR population to expected deaths in the general population (adjusted for age- and gender-specific rates) (Fig. 2). Thus, SMRs differ from hazard ratios as these compare the mortality of specific dcSSc patients with the mortality of other dcSSc patients. The general population mortality data were extracted from the 2014 Italian life tables  as Italian was the most prevalent nationality (22%) in the analyzed cohort from the EUSTAR database at the time of analysis . This was used as a mortality multiplier, sensitive to the fact that a patient is in a specific HAQ-DI state or any advanced organ state.
After calibration, a microsimulation structure was constructed, which enabled patients to move between the five HAQ-DI states, to develop major advanced organ involvement, and/or to enter a death state. An illustrative microsimulation was proposed for a patient to transition between the five HAQ-DI states, to develop lung/no lung involvement, and to enter a death state. The structure of the model is shown in Fig. 3. In that structure, a female patient starts in a specific HAQ-DI and lung state (defined by baseline characteristics). At the start of each cycle, based on her HAQ-DI state and lung state in the previous cycle, she could evolve towards another HAQ-DI and/or lung state, or die. The equations used to model the transitions in that illustrative model were: