The impact of tocilizumab on respiratory support states transition and clinical outcomes in COVID-19 patients. A Markov model multi-state study

Background The benefit of tocilizumab on mortality and time to recovery in people with severe COVID pneumonia may depend on appropriate timing. The objective was to estimate the impact of tocilizumab administration on switching respiratory support states, mortality and time to recovery. Methods In an observational study, a continuous-time Markov multi-state model was used to describe the sequence of respiratory support states including: no respiratory support (NRS), oxygen therapy (OT), non-invasive ventilation (NIV) or invasive mechanical ventilation (IMV), OT in recovery, NRS in recovery. Results Two hundred seventy-one consecutive adult patients were included in the analyses contributing to 695 transitions across states. The prevalence of patients in each respiratory support state was estimated with stack probability plots, comparing people treated with and without tocilizumab since the beginning of the OT state. A positive effect of tocilizumab on the probability of moving from the invasive and non-invasive mechanical NIV/IMV state to the OT in recovery state (HR = 2.6, 95% CI = 1.2–5.2) was observed. Furthermore, a reduced risk of death was observed in patients in NIV/IMV (HR = 0.3, 95% CI = 0.1–0.7) or in OT (HR = 0.1, 95% CI = 0.0–0.8) treated with tocilizumab. Conclusion To conclude, we were able to show the positive impact of tocilizumab used in different disease stages depicted by respiratory support states. The use of the multi-state Markov model allowed to harmonize the heterogeneous mortality and recovery endpoints and summarize results with stack probability plots. This approach could inform randomized clinical trials regarding tocilizumab, support disease management and hospital decision making.


Conclusion
To conclude, we were able to show the positive impact of tocilizumab used in different disease stages depicted by respiratory support states. The use of the multi-state Markov model allowed to harmonize the heterogeneous mortality and recovery endpoints and summarize results with stack probability plots. This approach could inform randomized clinical trials regarding tocilizumab, support disease management and hospital decision making.

Background
COVID-19 pandemic has brought major challenges for health care systems even in well-organized and resource-rich settings. Although the evidence for effective treatments are still missing, the number of observational and randomized trials on antiviral drugs and immune-active agents is rapidly increasing [1].
Several randomized controlled trials assessing antiviral drugs in COVID-19 patients have recently used an outcome that was suggested by the WHO R&D Blueprint expert group [2][3][4][5]. This outcome is a time to event categorical variable that considers the occurrence of death, recovery, but also the switch of oxygen support states, which represent well the patients' pathway in hospitals and provide information on their prognosis [6].
This outcome depicts progression and recovery of COVID-19 patients by their need of different respiratory support aids. A patient is classified according the following categories: 1. not hospitalized with resumption of normal activities; 2. not hospitalized, but unable to resume normal activities; 3. hospitalized, not requiring supplemental oxygen; 4. hospitalized, requiring supplemental oxygen; 5. hospitalized, requiring nasal high-flow oxygen therapy, noninvasive mechanical ventilation, or both; 6. hospitalized, requiring ECMO, invasive mechanical ventilation, or both; and 7. death [6].
The WHO R&D Blueprint score has never been used to assess the benefit of tocilizumab in the treatment of COVID-19 pneumonia. Tocilizumab is an extensively used immune active agent which showed discordant results in the reduction of IMV and mortality in hospitalized patients with severe COVID-19 pneumonia [7][8][9]. TESEO cohort showed a significant reduction in risk of death (aHR = 0.38, 95% CI:0.17-0.83, p = 0.02) and/or mechanical ventilation (aHR = 0.61, 95% CI:0.40-0.92; p = 0.02) in patients who received tocilizumab over the standard of care [10]. On the contrary, initial randomized clinical trials, did not show any benefit [7,11], while recent large randomized trials suggest reduced likelihood of progression to IMV and death in patients with severe COVID-19 pneumonia treated with tocilizumab [12][13][14]. In particular, an open-label randomized multicenter study conducted in Italy including 126 patients with PaO 2 /FiO 2 between 200 and 300 mmHg did not show the benefit on disease progression [7]. On the other hand, REMAP-CAP trial showed improved 90-day survival in critically ill patients receiving organ support who were treated with tocilizumab or sarilumab (HR = 1.61, 95% CI: 1.25-2.08) [12]. Moreover, RECOVERY trial reports that patients receiving tocilizumab had improved survival and were less likely to reach the composite endpoint of IMV or death (33% vs. 38%; risk ratio 0.85; 95% CI 0.78-0.93; p = 0.0005) [14].
These conflicting results could be justified by diverse inclusion criteria, small sample size, competing effect of other treatments such as glucocorticoids, but also by different stages of the disease when tocilizumab is administered [15].
The latter could be of paramount importance. Therefore, we hypothesized that a multistate Markov model applied in accordance to WHO R&D Blueprint indication, may help to better depict the impact of tocilizumab at any disease stage as captured by respiratory aid change during hospitalization. We hereby present a secondary analysis of a subset of hospitalized patients treated with tocilizumab previously described in TESEO cohort.
The aim of this observational study was to estimate the impact of tocilizumab administration on switching respiratory support states, mortality and time to recovery in hospitalized patients with SARS-CoV-2 pneumonia.

Population
This study included consecutive adult patients (�18 years) admitted to Infectious Disease Clinic of the University Hospital of Modena, Italy, from 21 February to 10 April 2020, with radiological findings suggestive for SARS-CoV-2 pneumonia and confirmed by PCR method on nasopharyngeal swab, belonging to the previously described TESEO cohort [10], but restricting the analyses to the Modena's cohort subset only. Inclusion and exclusion criteria were the same of the parental cohort [10]. Patients that did not satisfy criteria for severe pneumonia were excluded from the analysis. Each patient was followed up until discharge or death or up to 10 April 2020.

Definition of states characterising hospital clinical pathways
Patients received oxygen supply to target SaO 2 >90%. Following the approach suggested by WHO R&D Blueprint expert group, we defined five possible consecutive respiratory support states, that are hereinafter reported and referred as intermediate states: • State 1: no respiratory support (NRS); • State 2: oxygen therapy (OT); • State 3: non-invasive ventilation (NIV) or invasive mechanical ventilation (IMV); • State 4: OT in recovery; • State 5: NRS in recovery.
Two respiratory support states, such as non-invasive ventilation (NIV) or invasive mechanical ventilation (IMV) states were grouped to have a properly sized single category. Additionally, as per inclusion criteria only hospitalized patients were enrolled, we were not able to use the complete scale proposed by WHO R&D Blueprint expert group, which also comprises non-hospitalized patients.
Patients were admitted to the hospital at any state between 1 and 3 and during their hospital stay they did not necessarily switch through all states. The forward changes, which happen from any state between 1 to 3, indicated that patients were worsening, whereas those occurring from state 3 to 5 suggested that patients were improving.
Time was measured in minutes and reported in days, and was calculated as the elapsed time between the time the patient entered in one state and the time in which he/she moved forward to another one. Patients could leave the pathway for 2 possible reasons: 1. discharge, which can be at home, or at post-acute care facility; 2. in-hospital death.
Censoring was defined as the condition in which the patient was either still hospitalized on 10 April 2020 at 12:00 AM or transferred to another acute care facility. The representation of this sequence of states is reported in Fig 1. The main outcomes of the study were: progression from OT or NIV/IMV states to death, progression from OT to NIV/IMV, recovery from NIV/IMV state to the OT state. All the other changes in respiratory states were considered as secondary outcomes.

Pharmacological interventions for COVID-19 pneumonia
All patients were treated with standard of care in agreement with the Regional Guidelines of Emilia Romagna [16] regarding the treatment of COVID-19 that were continuously updated during the period of the study. The treatments consisted of: • Hydroxychloroquine (400 mg BID on day 1 followed by 200 mg BID on days 2 to 5 eventually adjusted for creatinine clearance estimated by a CKD algorithm); • Azithromycin (500 mg QD for 5 days) at physician's discretion when suspecting a bacterial respiratory superinfection; • Low molecular weight heparin for prophylaxis of deep vein thrombosis according to body weight and renal function.
• Lopinavir/ritonavir (400/100 mg BID) or darunavir/cobicistat (800/150 mg QD) for 14 days were used up to 18 March, when a clinical trial on the former did not show any benefit of protease inhibitors against the standard of care [3]. Other antiviral agents showing promising results such as remdesivir [17], was not used in our patients, due to its unavailability on the market.
Moreover, some patients were treated with immune-active agents, such as tocilizumab and glucocorticoids, that were not standard of care during observational period. Tocilizumab was offered by a treating physician to patients who experienced arterial oxygen saturation SaO 2 <93% and a PaO 2 /FiO 2 <300 mmHg in room air or a decrease in PaO 2 /FiO 2 > 30% in the previous 24 hours during hospitalization. All tocilizumab-treated patients provided a consent to participate in the study. Tocilizumab was administered intravenously (IV) or subcutaneously (SC) depending on hospital availability of the formulation at time of treatment [10]. IV tocilizumab was administered at the dose of 8 mg/kg of body weight (up to a maximum dosage of 800 mg) repeated twice, 12 hours apart, while SC formulation at the dose of 162 mg administered twice simultaneously, one per each thigh [18]. The patients who refused tocilizumab were treated with standard of care only.
Glucocorticoids were administered alone when tocilizumab was not available. On the contrary, patients underwent glucocorticoids after tocilizumab in the case of treatment failure, defined as lack of improvement of PaO 2 /FiO 2 after three days of tocilizumab administration, and these patients were included in the tocilizumab group. Patients who received anakinra or other immunomodulatory agents other than tocilizumab or glucocorticoids were excluded from the analysis.
All patients were evaluated possible pharmacologic drug-drug interactions before administration of any treatment and in particular cases, dosages were adjusted according to other comorbidities and co-medications. Patients were also carefully monitored for potential toxicities of each pharmacological agent.

Covariates
All the data was obtained from routinely collected medical records. The patients' full medical history, demographic and epidemiological data as well as the value of PaO 2 /FiO 2 at baseline were obtained at the hospital admission. The risk of multiorgan failure and mortality was assessed with standardized Subsequent Organ Failure Assessment (SOFA) score. Other covariates including age, gender, days from epidemic onset (defined as 21 February as the date of first COVID-19 case in the Modena province) were considered. Clinical data with patients' signs, symptoms, blood count, coagulation, inflammatory and biochemical markers were routinely collected and reported in an electronic patient charts. Other pharmacological interventions were glucocorticoids and lopinavir/ritonavir or darunavir/cobicistat. Use of glucocorticoids was considered as a potential confounder in adjusted models.

Statistical analysis
Continuous variables were described as mean or median and inter-quartile range (IQR), whereas categorical variables summarized as absolute and percentage frequencies. Comparison among patients receiving different types of immune-active agents were assessed by using the Yates chi-square test for categorical variables or the Wilcoxon test for quantitative variables.
Statistical analyses were carried out by using a continuous-time Markov multi-state model, which considered the previously defined sequence of respiratory support states. The transition probabilities, that is the probability to switch from a state to another one, were modelled according to an exponential distribution for time-to-event data, which takes into account censored follow-up times. Time was measured in days, considering the elapsed time in each transition across states. The Markov multi-state approach, compared to classical statistical methods for time-to-event data, has the advantage of allowing the joint analysis of length of care and incidence of clinical outcomes, such as the need for respiratory support or mortality, providing a complete assessment of COVID-19 inpatients' progress [19]. Under a statistical perspective this approach is also preferable as it can accommodate for competing outcomes, such as death or discharge, as well as for time-dependent multiple events, such as the need for respiratory support states [20].
At the beginning an unadjusted multi-state model analysis was carried out to assess the impact of tocilizumab on changes of state. Stratified predicted probabilities for all changes of states were then estimated, as well as the predicted lengths of stay in each state. These predictions were calculated for patients receiving tocilizumab since the beginning of the OT state and for patients never receiving it. A multivariable multi-state model was then considered, where the associations between the administration of tocilizumab and the changes of states were adjusted for the following variables: age, gender, SOFA score at hospital admission, time from the onset of COVID-19 in Modena and administrations of glucocorticoids and lopinavir/ritonavir or darunavir/cobicistat. Unadjusted and adjusted hazard ratios (HR) were used to measure the associations between tocilizumab administration and changes of state.
Treatment with azithromycin, hydroxychloroquine and low molecular weight heparin were not included in the multivariable model, as they were administered to all the patients. Time from COVID-19 onset in Modena was measured as the difference in days between the date of hospital admission and the 21 February 2020. The effect of age was reported as 5-years linear increment, the impact of SOFA score was described as 1-point linear increment, whereas the one related to time from COVID-19 onset as 1-week linear increment. The associations between the independent variables and the change of respiratory support states were assessed only for those transitions observed in at least 10 patients, whereas for the remaining transitions no independent variables were assessed in the models.
To properly assess the effect of tocilizumab on outcomes, its administration was modelled as a time-dependent covariate. This allowed to assess the effect starting from the exact time of first administration in each patient and to prevent the occurrence of immortal time bias. The effect of drugs administration was only estimated from the OT state onwards. Age, gender, SOFA score, time from the onset of COVID-19 in Modena and other pharmacological interventions were also evaluated as potential effect modifiers by adding interaction terms in the multivariable multi-state model.
To visualize the results simultaneously over time in a single informative figure, we used stacked probability plots [19]. These were calculated by using the unadjusted multi-state model's predicted prevalence of patients over time in each state, including discharge and mortality. The plots were referred to: patients who were treated with tocilizumab from the OT state onwards; patients who were not treated with tocilizumab during hospitalization.
HRs and predicted probabilities of change of respiratory support states were reported as point estimates and as 95% confidence intervals (95% CI). The number needed to treat was also calculated for those changes of state that describe a clinical worsening and for which a statistically significant results was observed. The goodness of fit of the models was assessed by comparing the observed and expected prevalence across time, for all respiratory support states, including hospital discharge and death. The analyses were performed by using R 3.6.3 statistical software (The R Foundation for Statistical Computing, Wien). The significance level was set to 0.05, except for HRs, which were tested at both p <0.05 and p <0.01 levels.

IRB approval and data information
Patients who received tocilizumab provided verbal, not written, informed consent because of isolation precautions. This procedure was approved by the Regional Ethical Committee of Emilia Romagna (protocol number: AOU 0018046/20). The same body also approved the retrospective analysis of the data and all patients provided informed consent to participate in the study. Oral consent was documented in the electronic patients' charts. Use of tocilizumab was based on clinical judgment and criteria reported in the TESEO study [10]. All patients with severe pneumonia admitted to hospital were considered for the analysis. Patients who consented to use of tocilizumab, but rejected to participate in the study were excluded from the analysis. Patients were admitted from 21 February 2020 to 10 April 2020. The retrospective data were fully anonymized and the only sensitive data was year of birth.

Transitions between states
Two hundred seventy-one patients were included in the analyses. The majority of patients were males (75.6%), the median age was 63 years. The overall follow-up time was 3159.8 total person-days (average 11.7, range 0.
No difference was found in age and gender between groups. The demographic and clinical characteristics of patients are reported in Table 1. Based on the observed data, some differences could be noticed in the two treatment groups.
Amongst patients who received tocilizumab 9 (7.4%) have died, whereas 38 (25.5%) patients have died without receiving tocilizumab. Considering patients who were given tocilizumab during NRS or OT states, 31 (32.3%) reached the NIV/IMV state. Conversely, amongst those who started their hospitalization in the NRS or OT states and were not given tocilizumab in these states, 69 (40.1%) underwent NIV or IMV.
The prevalence of patients in each respiratory support state was estimated in the stack probability plot, comparing people treated with and without tocilizumab since the beginning of the OT state (Fig 2).

Adjusted effect of tocilizumab
The results of the multivariable multi-state model, adjusted for age, gender, SOFA score, time from the onset of COVID-19 epidemic and administrations of glucocorticoids and lopinavir/ ritonavir or darunavir/cobicistat, were reported in Table 3. This analysis was carried out considering 260 patients with complete data. As in the unadjusted analysis, there was a positive effect of tocilizumab on the probability of moving from the NIV / IMV state to the OT in recovery state (HR = 2.6, 95% CI = 1.2-5.2). Furthermore, a reduced risk of death was observed in patients in NIV / IMV state (HR = 0.3, 95% CI = 0.1-0.7) or in OT (HR = 0.1, 95% CI = 0.0-0.8) treated with tocilizumab.
Time from the onset of COVID-19 epidemic was not associated with any risk of change of respiratory support state, whereas higher age, male gender and higher SOFA score were Interval. Predicted probabilities and HRs were estimated by using a Markov multi-state model. The predicted probabilities were referred to: patients who were treated with tocilizumab from the OT state onwards; patients who were not treated with tocilizumab during hospitalization. This analysis was not adjusted for confounders. # = statistically significant at 95% confidence level (p < 0.05) ## = statistically significant at 99% confidence level (p < 0.01).
https://doi.org/10.1371/journal.pone.0251378.t002 The predicted prevalence of patients in each respiratory support state was estimated in the stack probability plot, comparing patients who were treated with tocilizumab from the OT state onwards (left plot) and patients who were not treated with tocilizumab during hospitalization (right plot). Color code: black = death; blue = discharge; red = non-invasive ventilation/invasive mechanical ventilation; orange = oxygen therapy; green = no respiratory support. The strikethrough area refers to NRS or OT, whereas the "clean" areas refer to NRS or OT in recovery.
associated with several changes of respiratory state, all of them in the direction of a worse prognostic pathway (S1 Table). We also observed an increased mortality in patients treated with glucocorticoids when in the OT state and in patients treated with lopinavir/ritonavir or darunavir/cobicistat when in the NIV / IMV state (S1 Table). There was no evidence of modification of the effect of tocilizumab, due to age, gender, SOFA score, time from Covid-19 onset and other pharmacological interventions, based on significance testing of interaction terms.

Discussion
This prospective observational study showed the beneficial effect of tocilizumab in reducing the risk of death in patients with COVID-19 pneumonia. The reduced risk of death was larger at the initial oxygen support state, although this finding was based on a smaller number of patients, and persisted throughout advanced NIV/IMV support, respectively of 88% when used in OT and 70% in the NIV/IMV. Furthermore, patients receiving tocilizumab had an increased probability of recovery when in NIV/IMV state.
The reduced risk of death was previously described by our group in the context of the larger TESEO cohort [10] and we confirm, with a more detailed statistical approach, the lower mortality rate.
A systematic review and meta-analysis of 10 observational studies including 1358 patients demonstrated that mortality was 12% lower for COVID-19 patients treated with tocilizumab compared to patients who were not treated with tocilizumab. The number needed to treat was 11, suggesting that for every 11 COVID-19 patients treated with tocilizumab 1 death was prevented.
Moreover, encouraging results from large randomized clinical trials [12,14] could provide a greater level of evidence and recommendations in the current guidelines for the treatment of COVID-19 pneumonia.
This multi-state model allowed to depict a different beneficial effect size of tocilizumab according to different respiratory support states. We observed that tocilizumab may have had

PLOS ONE
The impact of tocilizumab on respiratory support states transition in Markov models a higher impact when used in the early disease stage, not only during IMV and pronounced cytokine storm. Tocilizumab blocks IL-6 receptors and along with TNF-alpha and IL-1 is involved in the pathways of cytokines' cascade [21]. In the late phase of acute respiratory distress syndrome (ARDS), when also other immunological factors are involved in pathogenesis of severe COVID-19 pneumonia, this agent might not be sufficient to face an uncontrolled inflammatory response [22]. The multivariable model did not include hydroxychloroquine, azithromycin and low molecular weight heparin because used as standard of care in all patients. In consideration of the fact that lopinavir/ritonavir or darunavir/cobicistat were administered to all the patients only before 18 March 2020, it was also decided to correct the multivariate analyses according to date of onset of the epidemic. Given the design of the study, we could not assess the role of tocilizumab without standard of care, therefore, we argue that tocilizumab alone may not be sufficient as a treatment of severe COVID-19 pneumonia. It is out of the purpose of this study to describe adverse events potentially associated with tocilizumab use, as being previously described in detail in the TESEO cohort [10].
The major confounders analyses in this study played a significant role in the progression of respiratory support states and therefore prevent from the most relevant sources of confounding bias. As previously described, age was a major driver of mortality in any disease stage [23]. Interestingly, age also had a negative impact in disease recovery. This is clinically evident in the loss of functional ability related to the intense catabolic damage and lean mass loss with progression of frailty [24]. Various studies reported a higher mortality rate among males [23,25]. However, the underlying mechanism is still unknown. It has been suggested that steroids and activity of X-linked genes can modulate the immune response to viral infection [26,27]. Additionally, ACE2 receptors might have a higher expression in men [28].
SOFA was included in the model as an indicator of hypertension, cardiovascular disease, chronic renal insufficiency or cancer. These disease conditions are potentially associated with organ failure during the acute inflammatory states induced by COVID-19 and, as expected they interfered with both intervention and outcome risks.
As previously described by our group, patients' trajectories across subsequent respiratory support states well depict the natural history and prognosis in hospitalized patients with COVID-19 pneumonia [submitted paper]. We propose a multi-state model that predicts a daily evolution of respiratory-support aids in patients with COVID-19 pneumonia and can be used to identify patients who will likely benefit most from immuno-active drugs. The application of this Markov model has been widely used in the setting of chronic disease conditions and have been rarely applied to the rapidly evolving clinical evolution of COVID-19 pneumonia [19].
Our statistical approach accommodates the primary analysis with a stacked probability plot of the major events such as being discharged alive and death. This model is a powerful tool to describe clinically opposite endpoints (mortality, discharge) as well as endpoints influencing hospital capacities (duration of hospitalization) simultaneously over time. In particular, the stack probability plots visualize results over time in a single informative plot. These figures allow the joint evaluation of multiple endpoints by visual comparison of depicted areas.
As described by Von Cube et al [6], the heterogeneity in choice of endpoints reflects the manifold potential treatment effects, the various different decision makers involved and the different patient populations under study. Nonetheless, harmonization of the different endpoints is an essential step to allow comparison of results from clinical trials and fasten decision making.
We acknowledge some limitations of our study. Firstly, we cannot exclude residual confounding bias by indication due to patients' characteristics that we did not adjust for, such as body mass index or comorbidities. However, by including clinically meaningful covariates such as age, gender and risk of multi-organ failure (SOFA) in the multivariable model, these bias can only be caused by other characteristics and should be of acceptable magnitude. Moreover, we also controlled for confounding related to time from COVID-19 onset, since in this study there were patients treated at the very beginning of the outbreak as well as several weeks later. Secondly, results may be partly influenced by the assumptions of the multi-state model, namely those concerning the onset and the duration of clinical effect of tocilizumab. Thirdly, this study was performed in one single tertiary centre in Italy and may thus not reflect general management of SARS CoV-2 pneumonia in other settings. Fourthly, NIV and IMV were grouped together in order to have properly sized single category, but this choice did not permit to distinguish the escalation from NIV to IMV and de-escalation from IMV to NIV as important steps in management of severe COVID-19 pneumonia. Fifthly, multistate models have the limitation that they do not account for patient triage possibly due to ICU congestion. Additional research is needed to avoid selection bias arising from this special clinical situation. This might be due to different epidemic environments over countries and due to heterogeneity in hospital health care capacity.
To conclude, we were able to show the positive impact of tocilizumab used in different disease stages depicted by respiratory support states. The use of the multi-state Markov model allowed to harmonize the heterogeneous mortality and recovery endpoints and summarize results with stack probability plots. This approach could inform randomized clinical trials regarding tocilizumab, support disease management and hospital decision making.
Supporting information S1