We built a microsimulation model of extreme preterm infants (gestational age less than 29 weeks) to estimate the long-term clinical and economic outcomes of patients stratified on gestational age at birth, BPD severity at 36 weeks, and highly associated complications. Our study combined the best available evidence on the BPD risk among extreme preterm infants, BPD severity-adjusted risk of developing major complications, quality of life and health system cost for each complication in order to estimate the lifetime burden associated with BPD.
We performed a targeted literature search without any restrictions on study design or a publication date. We included studies that were based on a single centre and those were based on the large national or regional neonatal databases, including (Canadian Neonatal Network (CNN), EPICure (UK), NICHD Neonatal Research Network Generic Database (USA), Vermont Oxford Network (USA), Kaiser Permanente Medical Care program (USA), Neonatal Research Network of Japan, Epipage 1 & 2 (France), AOK National Insurance entries (Germany). We used the following free-text terms: ‘bronchopulmonary dysplasia’ OR ‘BPD’ AND ‘preterm’ OR ‘premature birth’ OR ‘very low birth weight’ OR ‘extremely low birth weight’. Free-text terms were ‘long term’ AND ‘outcomes’ OR ‘follow-up’ OR ‘complications’. The bibliographies of included studies and pertinent reviews was also be hand searched for relevant studies.
There is limited population-based data on BPD status and associated complications according to even basic cohort stratifications such as gestational age or birth weight. The CNN collects nationwide real-world data and publically reports aggregated means on discharge status, mortality rates, and some complication outcomes according to patient characteristics that would be captured in a high-level chart review of the index admission (5). CNN data is limited in its tracking of patients post-discharge, but does provide a helpful starting point on the distribution of patients discharged from hospital. The published literature provides additional neonatal patient information on discharge rates, ventilation duration, long-term complication rate, and health care utilization (2, 6, 7). However much of the data were reported as aggregated means from small cohorts, and only one study stratifies according to BPD severity (2). The dearth of information on risk-adjusted outcomes and distributions around aggregated means necessitated a simulation approach to combine data sources and estimate patient distributions across age and risk categories (in our case, BPD status).
A microsimulation approach was used to facilitate the creation of sample distributions of patients at birth using the age-adjusted risk of death and the risk of developing BPD. BPD status was divided into severe, moderate, and mild. While there is debate within the field about the precise thresholds to classify severity of BPD, the relevant source data we utilized distinguished mild, moderate, and severe BPD according to a fraction of inspired oxygen (FiO2) of 0.21, < 0.30, and ≥ 0.30 at 36 week post-menstrual age respectively (2). We combined no BPD status with mild since there is minimal evidence suggesting a significant difference in patient outcomes between these two groups among extreme preterm infants. Available data provides mean estimates of mortality during admission and the distribution of BPD severity at 36 weeks of age and includes confidence intervals for chronic lung disease status and increased risk of death according to specific age thresholds (5). However, the reported evidence does not allow us to directly observe the mortality rates by both scales of BPD severity and gestational age categories necessary to distinguish outcomes at the granularity needed for decision modeling. We used the available confidence intervals to bind our analysis to avoid extreme estimates, but the model requires parameters categorized according to BPD severity and gestational age at birth the disease model to sufficiently address clinical and decision questions in the future.
We performed a first-order microsimulation and generated a sample distribution of patients with an assigned BPD severity status at 36 weeks from a ranked sample selection. We then applied a different weight to a fixed mortality rate according to severity group (severe, moderate, mild) that assumed infants with more severe BPD status are more likely to die prior to formal diagnosis at 36 weeks. The starting rates were sampled from published mortality rates (5), with the weights applied over 5,000 sample iterations, and then compared to the aggregated mortality rate obtained from the published literature. Weighted mortality rates that produced aggregated death rates that fell within the published confidence range were saved by the sampling algorithm, meaning as the sampling process was repeated the potential range of potential mortality rates within each BPD severity group became increasingly clustered on values that would produce aggregated rates that matched the real-world evidence. This process generated a range of potential risk-adjusted mortality rates for each BPD severity category. Using these risk-adjusted mortality rates, the ranked sampling process was repeated to draw patients by BPD severity status adjusted for mortality during admission that would skew the distribution of patients still alive at 36 weeks. (8) The result served as a prior patient distribution of BPD at birth and at 36 weeks adjusted for differential mortality risk; with the former having a disproportionately higher number of more severe cases that were more likely to not survive to discharge. The simulated patient distributions at 36 weeks were validated against the published distributions and 94% of our sampled mortality rates fell within the reported 95% confidence range (5). We followed the same sampling process for three different gestational age groups—less than 25 weeks at birth, 25–26 weeks, and 27–28 weeks—and combined the age groups into a final population based on the CNN reported age distributions.
Birth weight was excluded as risk factor due to minimal evidence on its impact on relative risk across both axes (mortality and BPD severity). Several sub-analyses were performed to test the impact of non-Gaussian distributions around the CNN data’s reported mean. Mortality risk by BPD status was found to significantly impact the final sample distribution. We concluded, on the basis of expert opinion and the published literature (4,6,7), that the mortality risk is non-normally distributed by patient characteristics and is very likely to be a competing risk of BPD status at 36 weeks. Consequentially, we applied this non-Gaussian distribution of mortality risk to the final model by applying BPD-correlated ranks to the random effects model to estimate mortality and back transformed to a uniform distribution to apply to the final joint density function. (9)
Using the first-order probability density function approach described above, a final patient population of 10,000 simulants were used to represent the distribution of gestational age at birth and BPD status. The population represents infants that survived to discharge. We modelled the lifetime trajectory of patients using an annual cycle. Post discharge, patients could develop major complications highly associated to preterm infants and differentially distributed according to BPD. We included the following major complications: respiratory illnesses requiring rehospitalization with variable length of stay, clinic visits for wheezing episodes, asthma exacerbations or psychiatric illnesses such as ADHD, developmental delay which can be global or specific (global developmental delay is defined in reference to infants and preschoolers, ages 0–5 years, who present with delays of 6 months or more, in two or more of the following developmental domains: gross/fine motor, speech/language, cognition, social/personal and daily living activities while specific developmental delay (SDD) refers to age inappropriate performance in a specific area), neurological impairment (cerebral palsy), hearing impairment, vision impairment, and chronic pulmonary hypertension. (10) We performed a targeted literature search to identify the risk of developing major complications according to BPD status and average age of diagnosis in order to derive time-dependent transition probabilities for each complication.
Costs and health utilities associated with each complication were applied if the infant acquired the complication, with all costs and utilities adjusted for half-cycles in the cycle they acquired the complication and the year of death. We estimate healthcare cost from a perspective of the Canada’s publicly funded healthcare system; only costs that are borne to the health system are incorporated. Cost information was taken from published Canadian studies that specified the cost burden attributed to BPD or a specific complication (see Table 2). Costs directly pertaining to BPD included initial admission, oxygen therapy, any direct long-term costs that were identified as treatment for BPD (2), and clinic visits and hospital admissions within the first two years. Our final cost estimate for burden of BPD included cost estimates for major complications linked to BPD. When applicable, costs were time-dependent to reflect changes in healthcare utilization of a patient’s lifetime. We estimated quality-adjusted life years (QALY) by multiplying the expected life expectancy and respective health utility values. There is evidence on patient health utilities for BPD patients that differentiate by severity (11) however adjusting utilities according to an underlying BPD status differentially by combinations of complications presented significant challenges without much existing evidence to validate this approach. Consequentially, we used two different health utility measures to account for differences in expected patient utility based on the combination of complications developed over their lifetime. The base case approach assigns patients a utility score according to their BPD severity (Mild being the relative utility reference) with the complication a patient develops that has the lowest relative utility estimate being applied as additional disutility. An alternative approach treated each complication acquired as compounding a patient’s disutility and were applied multiplicatively, thereby presenting a more extreme interpretation of quality of life. Costs and utilities were discounted at 1.5% in accordance with evaluation guidelines (12). Costs are presented in 2018 Canadian dollars.
Table 1
Input Parameters for Risk of Event/Complications*
| | BPD Severity Status |
| | Mild | Intermediate | Severe |
| | Mean | SE | Mean | SE | Mean | SE |
Event | 2-year Clinic Visits | 3.1 | 1.9 | 28.5 | 17.1 | 28.5 | 17.1 |
Hospital Admissions in Year 1 | 3.4 | 1.7 | 4.3 | 2.2 | 7.7 | 3.9 |
Hospital Admissions in Year 2 | 2.3 | 1.2 | 2.6 | 1.3 | 3.6 | 1.8 |
Length of Hospital Stay (Year 1 & 2), days | 18.8 | 2.3 | 8 | 2.2 | 10.9 | 2.0 |
Risk of Complication | Neuro-impairment | 34.4% | 6.9% | 58.4% | 11.7% | 78.9% | 15.8% |
ADHD** | 8.6% | 1.7% | 5.8% | 1.2% | 6.5% | 1.3% |
Neurodevelopmental Delay | 14.3% | 2.9% | 23.4% | 4.7% | 21.5% | 4.3% |
Asthma | 35.9% | 7.2% | 35.9% | 7.2% | 34.0% | 6.8% |
Hearing Impairment | 5.7% | 1.1% | 14.3% | 2.9% | 10.5% | 2.1% |
Low Vision | 9.1% | 1.8% | 20.8% | 4.2% | 21.1% | 4.2% |
Pulmonary Hypertension | 0.0% | 0.0% | 4.7% | 0.9% | 6.5% | 1.3% |
* Sources: 2,3; **ADHD = Attention Deficit Hyperactivity Disorder |
Table 2
Input Parameters Used for the Model
| | Proportion | Mean cost (C$, 2018) | SD | Source |
Costs |
Preterm Infant Index Admission | Extreme preterm (< 28 weeks) | | 59,508 | 10,473 | 12 |
Preterm Infant Long Term Costs | Extreme preterm (< 28 weeks) | | 12,884 | 2,150 | 2 |
BPD (By severity) -Index Admission Cost | All BPD | | 105,131 | 18,393 | 3 |
Mild | | 75,928 | 15,154 |
Moderate | | 113,892 | 19,504 |
Severe | | 141,148 | 27,823 |
BPD Home Care Costs | All BPD | | 475 | 81 |
BPD Annual Costs (Post Discharge) | All BPD | | 15,628 | 3,056 |
Annual management costs for general population (or preterm infant specific) with: |
ADHD | Drug Costs | 30% | 827 | 144 | 13 |
Psychological/Behavior therapy | 15% | 2,880 | 523 |
Combination (Medication + Counseling) | 32% | 3,707 | 612 |
No treatment | 23% | 189 | 34 |
Asthma | | | 3,925 | 709 | 14 |
Hearing impairment | Cost of hearing aid device | 76% | 1,549 | 295 | 15 |
| One time cost of cochlear implant | 24% | 46,555 | 8,376 |
| Post implant follow up year 1 | 4,411 | 825 |
| Post implant follow up year 2 | 1,766 | 285 |
| Post implant follow up year 3 | 1,255 | 216 |
Retinopathy | | | 3,914 | 750 | 16 |
Pulmonary Hypertension | Total healthcare costs | | 18,882 | 3,449 | 17 |
Neuro-impairment | Without technical assistance | 80% | 11,900 | 1,976 | 18 |
| With technical assistance | 20% | 46,603 | 8,791 |
Developmental Delay | Total healthcare cost | | 10,534 | 1,861 |
Utilities |
Preterm BPD | Mild BPD (reference category) | | 1.00 | | 19 |
Intermediate BPD | | 0.80 | 0.20 |
Severe BPD | | 0.50 | 0.10 |
ADHD | Average | | 0.71 | 0.25 | 20 |
Asthma | Average | | 0.89 | 0.09 | 21 |
Hearing impairment | Average across the severity (mild, moderate, severe hearing loss) | | 0.62 | 0.20 | 22 |
Retinopathy Of Prematurity | Bilateral threshold retinopathy of prematurity | | 0.60 | 0.10 | 23 |
Pulmonary Hypertension | Average | 0.71 | 0.14 | 24 |
Neuro-impairment | Cognitive impairment in preterm sample | 0.64 | 0.33 | 25 |
Developmental Delay | Average across levels | | 0.42 | 0.41 | 26 |
Note: BPD = Bronchopulmonary dysplasia; ADHD = Attention Deficit Hyperactivity Disorder |
Population age-dependent mortality risk was derived from Statistics Canada (13). Adjusted relative risk of mortality following discharge for extreme preterm infants was applied to all patients based on best available evidence (14). There is currently no long-term study that observes relative risk of mortality according to BPD status, so mortality risk is assumed to vary by gestational age at birth in this model.
All results are presented as the mean probabilistic output following a second-order Monte Carlo of 10,000 iterations wherein all parameter values were varied according to their distribution. Each second-order run utilized a starting distribution of patients distributed according to their BPD status at discharge. 1,000 first-order runs of the patient distribution at discharge were run for every second-order iteration. The model was developed using R (R Foundation for Statistical Computing, Vienna, Austria).