In silico clinical trials for pediatric orphan diseases

To date poor treatment options are available for patients with congenital pseudarthrosis of the tibia (CPT), a pediatric orphan disease. In this study we have performed an in silico clinical trial on 200 virtual subjects, generated from a previously established model of murine bone regeneration, to tackle the challenges associated with the small, pediatric patient population. Each virtual subject was simulated to receive no treatment and bone morphogenetic protein (BMP) treatment. We have shown that the degree of severity of CPT is significantly reduced with BMP treatment, although the effect is highly subject-specific. Using machine learning techniques we were also able to stratify the virtual subject population in adverse responders, non-responders, responders and asymptomatic. In summary, this study shows the potential of in silico medicine technologies as well as their implications for other orphan diseases.

Scientific RepoRTs | (2018) 8:2465 | DOI: 10.1038/s41598-018-20737-y parameter sets can be created in silico, a large cohort of virtual subjects can be established, overcoming the limitation of small patient cohorts in rare diseases. Furthermore, the challenges associated with (pediatric) randomized clinical trials are avoided by using the in silico model to generate a unique paired data set of both treated and non-treated virtual subjects.
In this study we aim to combine data-driven and mechanistic modeling approaches to illustrate how in silico models can contribute to unraveling the mechanisms underlying orphan diseases, with a particular focus on congenital pseudarthrosis of the tibia (CPT) associated with mutations in Neurofibromatosis Type 1 (NF1). CPT is a rare disease which is characterized by a progressive bowing of the tibia and develops into spontaneous fractures in the distal third of the tibia 12 . Typically the bone regeneration is insufficient, resulting in a pseudarthrosis, which is treated by physical excision of the abnormal bone tissue (called the 'hamartoma') and internal or external fixation 12 . In recent clinical practice, rhBMP-2 or rhBMP-7 (recombinant human bone morphogenetic protein) is used to enhance the surgical outcome. However, the effectiveness of BMP therapies is still highly debated 13,14 and the FDA has issued a warning against the use of BMP in skeletally immature patients 15 . As such, the goal of this study is to assess the effect of BMP treatment in an extensive, investigative in silico clinical trial using a previously developed mechanistic model of murine bone regeneration 16 and advanced machine-learning techniques to mine the large set of virtual subject data.

Results
The virtual clinical trial predicted the efficacy of treatment. The degree of severity of CPT was assessed by the complication index (CI), a linear combination of three phenotypic CPT symptoms (see Methods). A parameter combination yielding a small CI value reflects a fairly normal fracture-healing process with only minor aspects of CPT (if any). Conversely, a parameter combination resulting in a large CI value represents severely impaired fracture healing reminiscent of the CPT phenotype. The virtual clinical trial predicted a wide range of CI values, corresponding to different degrees of severity by which the pseudarthrosis presents itself, for both the treated and the non-treated group of subjects ( Fig. 2A). The group of treated subjects was characterized by a significantly lower mean CI value (Fig. 2B, *p < 0.01) and a higher amount of bony unions (Fig. 2C), indicating a beneficial effect of BMP treatment.
The virtual clinical trial stratified subjects. To see whether the whole virtual subject population responded similarly to the BMP treatment, we plotted for each subject the results of the treatment versus no treatment (Fig. 2D). Using an arbitrary cut-off of CI value of 0.5, the virtual subjects clustered into four distinct groups: adverse responders (a low CI when untreated and a high CI when treated, red); non-responders (a A previously developed mechanistic model of murine bone regeneration was adapted to describe the aberrant cellular behavior of cells affected by NF1 mutation by altering the corresponding parameter values (inner, blue circle: in silico model).
As an example, the inner circle shows how a conceptual model (1) of cartilage generation can be translated in a schematic (2) and mathematical representation (3). By altering the parameter values, various physiological and pathological behaviors can be represented (indicated in green). From the altered in silico model, a cohort of 200 "virtual subjects, " each representing a different parameter combination, was generated. For each virtual subject the bone-healing process was simulated both with BMP treatment and without BMP treatment (outer, green circle: in silico clinical trial).
Interestingly, unsupervised machine-learning methods (Ward hierarchical clustering) computationally confirmed the existence of three classes (Fig. 2E): non-responders (blue), asymptomatic (green), and responders (black). The algorithm did not define a separate group of six adverse responders as distinguished by the manual method. However, five of these subjects clustered together in the asymptomatic group, and one was placed among the non-responders group (blue) by the algorithm. The Ward hierarchical clustering algorithm also classified three responder subjects (black) among the asymptomatic (green), which is different from the manual classification. Note that these four subjects differently classified by the manual and computational methods all lie very close to the previously defined arbitrary cut-off (Fig. 2D).
The virtual clinical trial identified potential biomarkers. Next, we investigated whether the manually defined subject classes have distinctive characteristics that can be used as potential biomarkers to predict the efficacy of BMP treatment. We focused here on the difference between the responders (black) and the non-responders (blue), as these subjects can benefit most from BMP treatment. We first calculated the Spearman correlation matrix of the dataset including the 8 NF1 parameters as well as fibrous, cartilage, and bone tissue fractions at different time points to check for highly correlated (redundant) attributes (Fig. 3). Importantly, the results showed that the NF1 parameters were uncorrelated. As expected, the CI value was negatively correlated with the bone tissue fraction and positively correlated with the fibrous tissue fraction. Similarly, some NF1 parameters were correlated with certain tissue fraction outcomes (e.g., A f0 , which describes the rate of fibroblastic proliferation (see Methods), with the amount of fibrous tissue), and the tissue fractions were also not independent.
Using the 8 NF1-related parameter values in data-driven modeling, we determined which parameter was most important in predicting stratification between paired subject groups ( Table 1, ROC-curves can be found in the supplementary material). Note that not all comparisons are balanced, particularly for the adverse responders with only 6 subjects. Across all paired groups, P mc , Y 11 , and Y 3cb are important predictors (Table 1, Supplementary  Figure 1), with mean accuracy values >0.72. P mc represents the rate of cartilage formation, Y 11 signifies osteogenic differentiation, and Y 3cb denotes endochondral ossification (see Methods). Importantly, when similar analyses are performed on scrambled data, a mean accuracy of 0.48 was obtained, and no clear predictors were found (Supplementary Figure S1D-F). These data independently validated the accurate stratification of the subject groups.
In order to analyze the data on CI in more detail, we looked at the different components constituting the CI value, e.g., the amount of fibrous tissue after 49 days and the amount of non-unions (Fig. 4A). Interestingly, the box plots show that there was a significant difference in the amount of fibrous tissue present at day 49 after treatment for the different subject classes, with significantly less fibrous tissue present in the responding subjects. Moreover, the fraction of non-unions is also lower in the responding subjects (Fig. 4B).

Discussion
The clinical literature reports conflicting evidence about the effectiveness of BMP in treating CPT 13,14,17,18 . It is difficult to compare these clinical studies and draw solid conclusions for various reasons. First the surgical procedures and fixation methods are very heterogeneous, varying between Ilizarov fixation 17 , intramedullary fixation 14 , autologous grafting 18 , Masquelet technique 13 , or a combination thereof. Second, some clinical studies are performed with rhBMP-7 18 whereas others use rhBMP-2 14 . Third, the delivery (e.g., mixed with collagen and allograft bone 18 or the use of a collagen sponge alone 14 ) and the dosage of the rhBMP (e.g., from 3.5 18 to 12 mg 14 ) varied between studies. Fourth, with CPT being an orphan disease, the few case studies with limited follow-up of refracture result in a small amount of patient-specific data. Finally, and most importantly, the treated CPT patients each have a specific case history with previous treatments, including multiple surgeries and rounds of BMP treatment, which makes conclusions hazardous.
Since virtual clinical trials alleviate some of the above mentioned problems, we combined in this study mechanistic modeling with data-driven modeling in an investigative in silico clinical trial to determine the (beneficial)  Table 1. The NF1-related parameters with greatest predictive power to differentiate the paired subject groups in the virtual clinical trial a . a Predictors are as defined in the Methods and based on mean arbitrary units >50 with a random forest, binary classification model. Mean accuracy ± standard deviation is given in parentheses and were measured across 100 trained random forest models.
effect of BMP treatment on fracture healing in NF1 subjects. A large in silico trial was designed, simulated and analyzed, based on a previously established computational model of murine bone regeneration. Notably, the computational model has been previously corroborated in various pathological conditions including small bone defects 19 , large bone defects 20,21 and CPT in NF1 subjects 16 . Although the parameters and the in silico clinical results are related through mathematical equations, the set of partial differential equations is very extensive and has many non-linear interactions between the variables. As such it is impossible to intuitively predict how the parameter values will influence the in silico predictions, here the bone healing outcome. Therefore, in this study, we decided to generate an extensive data set in silico and mine it using machine learning approaches. An important simulation outcome is the CI value, which reflects the degree of severity of CPT. We found that the CI value was significantly lower with BMP, compared to no treatment, for the 200 subjects in our in silico clinical trial (Fig. 2). Our computational model predicted a wide range (0.01-0.89) of CI values for both treated and untreated conditions, reminiscent of the substantial variation in severity of pseudarthrosis, which suggests that not all subjects benefit from treatment. The large range of predicted CI values for the non-treated subjects is remarkable since the parameter values of the eight factors were mostly varied in one direction with respect to the baseline value, thereby biasing the design to a CPT phenotype. By manual thresholding the CI values, four distinct subject groups were identified: adverse responders, responders, non-responders and asymptomatic. Unsupervised Ward hierarchical clustering confirmed the existence of all classes except the adverse responders, although these subjects also clustered together. The existence of four distinct subject groups may explain the heterogeneity in the reported clinical outcomes mentioned above, although further research is necessary to confirm this subject stratification and find biomarkers thereof, in particular for the small subpopulation of subjects predicted to negatively react to the BMP treatment.
It is also noteworthy that the in silico responders (n = 47) have a wide range in the amount of fibrous tissue at day 49 (Fig. 4A). The presence of fibrous tissue in the regenerating tissue may correlate to inferior mechanical properties and a higher risk for refracture. As such, the computational results seem to confirm the clinical observation that BMP induces healing but with a substantial refracture risk although this biomechanical hypothesis should be investigated further.
In addition to subject stratification, the data-driven modeling techniques allowed us to identify three NF1 parameters that can accurately predict stratification among paired subject groups (Table 1). With further validation, these parameters might form the basis of biomarkers for successful BMP treatment. The rate of cartilage formation (P mc ) is a distinguishing factor between the non-responders and the responders. This result corresponds to the findings of a previous computational study 16 where it was shown that P mc is a critical determinant for the CI value due its important role in the endochondral ossification pathway. Interestingly, Y 3cb , a factor representing endochondral ossification, was also identified as a predictive NF1 parameter. It should be noted that P mc and Y 3cb are not correlated.
While this study focused specifically on the design, execution and analysis of an investigative in silico clinical trial for CPT, a pediatric orphan disease, others have used in silico clinical trials to investigate different clinical  23 . Besides refining and partially replacing animals or humans in a clinical trial, in silico models can also generate virtual patients with particular characteristics to augment clinical trials, which is particularly relevant for trials in children or rare diseases 24 .
In addition to their value in studying clinical problems, in silico clinical trials offer huge potential in the design-discovery phase of biomedical product development. These virtual studies can identify biomarkers, optimize the dosage and duration of therapeutic interventions and patient stratification as well as perform risk and safety assessments required for regulatory applications 10 27 . Importantly, the in silico metabolic model has been shown to adequately represent glucose fluctuations in type 1 diabetes and has been accepted by the Food and Drug Administration as a substitute to animal trials in the preclinical testing of closed-loop control algorithms. Notably, because of the increasing interest in the use of in silico models, frameworks and guidelines are being developed to thoroughly investigate the reliability of in silico model predictions for regulatory applications 28 and in silico clinical trials 6,29 .
Although the in silico clinical trial presented here offers much insight into treatment of CPT, we recognize several limitations. Firstly, the in silico clinical trial considers a cohort of (only) 200 murine subjects in order to balance between statistical accuracy and computational feasibility (one simulation requires on average 48 hours computing time on one node of a high performance cluster). The virtual CPT subjects were based on a murine model of bone regeneration, including a generic murine-based geometry, due to the lack of experimental human data and to reduce computational demands. In silico simulations have shown, however, that increasing the fracture geometry to human dimensions will proportionally increase the fracture healing time (unpublished data). Secondly, the virtual subjects were uniformly sampled across the parameter space which resulted in a large range of phenotypes and inherent variability in the data. From the additional data-driven analyses, four subject groups were identified but these were not similar in size. Moreover, the group size of the adverse responders was too small to draw any statistically significant conclusions. In the future, advanced algorithms, such as bootstrap aggregating, could be used to tackle the issues related to unequally sized clusters. A larger cohort of subjects could also be simulated, with an improved sampling of the group of adverse responders as well as a more refined definition of the manually identified groups. Thirdly, the degree of severity of CPT was assessed by the CI, previously defined by three phenotypic criteria: the absence of cortical bridging of the defect, the amount of fibrous tissue and the amount of fibroblasts 16 . The results of this clinical trial should be interpreted in light of this definition, which puts a lot of importance on the presence of fibroblasts and fibrous tissue.
Nevertheless, the combined mechanistic and data-driven modeling approach allows us to successfully perform and analyze an investigative in silico clinical trial, including the identification of subject subpopulations. An obvious next step would be to extensively characterize the real patient-specific parameter distributions and the corresponding distribution of clinical CI values and compare those with the simulated cohort. However, due to the particular nature of the disease (orphan, pediatric), it will be extremely difficult to obtain these distributions for a sufficiently high number of patients. As an intermediate proxy, we intend to characterize a number of in vitro properties such as trilineage differentiation capacity of NF1-haplodeficient cells and NF1-null cells from various donors. Another important step towards a clinical translation of these results is the in vivo validation of the in silico predictions in order to find out whether these subject subpopulations truly exist. Moreover, in case of discrepancies between the predicted and measured in vivo outcomes, the additional in vivo experimentation would allow for the identification of any shortcomings of the computational model, such as particular model simplifications (see ref. 16 for a more elaborate discussion on model simplifications).
In summary, in silico models and virtual clinical trials are increasingly being explored to tackle the challenges associated with pediatric diseases and/or orphan diseases. This study demonstrates how mechanistic and data-driven modeling are useful tools to simulate and mine data from in silico clinical trials and stratify subject populations to improve current treatment strategies. This robust modeling approach can also be employed for biomarker identification, optimization of the dosage or of the duration of the proposed intervention. Moreover, this study shows the vast potential of in silico medicine technologies as well as their implications for other orphan diseases.

Materials and Methods
Design of the in silico clinical trial. The set of 200 virtual subjects was generated from a previously established multiscale computational model of murine bone regeneration, which has been corroborated for small bone defects 19 , large bone defects 20 , and CPT in NF1 subjects 16 . Although the exact etiology of CPT is still unknown, 40-80% of CPT patients are carriers of a mutation in the NF1 gene. Also, it has been postulated that double inactivation of the NF1 gene is necessary for the development of pseudarthrosis 30,31 , and that the relative proportion of NF1-haplodeficient cells and NF1-null cells in the pseudarthrosis region will determine the healing outcome. To examine the effect of the NF1 mutation on bone-fracture healing, the parameter values of the factors describing the aberrant cellular behavior of NF1-haplodeficient and NF1-null cells were varied across a wide range in the computational model (Table 2)  lesion cells (c f,BC ) 12 , increased fibroblastic proliferation (A f0 ) 12 , increased fibroblastic differentiation (F 4 ) 32 , reduced osteogenic differentiation (Y 11 ) 12,33 , reduced endochondral ossification (Y 3,cb ) 12 , reduced cartilage formation (P mc ) 12 , increased fibrous tissue formation (P mf ) 32 , and increased angiogenic growth factor production (G gvc ) 12,34 . Importantly, the computational model focuses on the cellular and biochemical aspects of CPT and implicitly assumes that the fracture is sufficiently stabilized to allow bone formation. Hence, mechanoregulatory stimuli were not considered in the computational model.
The in silico clinical trial consisted of 200 virtual subjects for which the healing process was simulated both with BMP treatment and without treatment. The 200 virtual subjects were generated using the "Design of Experiments" (DOE) tool of JMP (SAS Institute Inc., Cary, North Carolina, USA) 35,36 . DOE (or experimental design) is a statistical method that generates an array of combinations of different parameter values within a predefined parameter space. A uniform design was chosen since the computational model of bone regeneration is highly non-linear. Uniform designs also cope well with the addition and removal of runs. The parameter space (Table 2, NF1 range) was determined by varying the parameter values of the eight factors mostly in one direction with respect to the normal case, in this way biasing the DOE design to a CPT phenotype. As such, each virtual CPT subject had a unique combination of eight parameter values, representing their intrinsic cellular and biochemical properties. Please note that the resection of the fibrous hamartoma, consisting of NF1-haploinsufficient and NF1-null cells, is not captured in the in silico model, or in other words, it models the worst-case scenario where NF1-null cells are still present in the fracture callus. The computer simulations were performed with these parameter combinations, utilizing high-performance computational resources provided by the VSC (Flemish Supercomputer Center).
BMP treatment was modeled using one generic osteochondrogenic growth factor (g bc ), which represents the effect of multiple growth factors present in the fracture callus and released from the BMP sponge as a clinical treatment. The influence of this generic growth factor on differentiation is either chondrogenic or osteogenic depending on the local oxygen tension. To model the gradual release of BMP from a BMP sponge, implanted directly after the occurrence of the fracture, a time-dependent, periosteal boundary condition was simulated (Eq. 1): is equal to the standard value of the osteochondrogenic boundary condition (see supplementary material).
The degree of severity of CPT is assessed by the complication index (CI, γ 7 ), which was previously defined in ref. 16 as a linear combination of the three most prevalent symptoms of CPT (a non-union (γ 6 ), harmatoma (γ 4 ), and presence of fibroblasts (γ 5 )) (Eq. 2): 3 (2) 7 4 5 6 γ γ γ γ = + + A parameter combination yielding a small CI value reflects a small degree of severity of CPT, or in other words, the fracture healing proceeds fairly normally. Conversely, a parameter combination resulting in a large CI value represents severely impaired fracture healing and is reminiscent of the CPT phenotype.

Statistical analysis.
For the statistical analysis of the CI values with and without BMP treatment, a two-sided, paired t-test was performed. To classify subjects based on their CI values with and without BMP treatment, we used the Ward hierarchical clustering algorithm (ward.D2, Euclidean distance) from the "stats" package implemented in R (version 3.3.1) 39 . To determine the prediction accuracy of the classification and the most important predictors, we created a machine-learning model with the "caret" 31 package employing the random classification algorithm from the "randomForest" package in R 40 . The random forest model was trained on 75% of the data with 10-fold cross-validation 41 . The accuracy was assessed on the remaining 25% of the data. To get sufficient precision, the random forest model was trained 100 times with random splitting of data in training and testing data set every time. The same analyses were also performed on scrambled data, obtained with the randomization algorithm from the "picante" package 42 . The script used to perform the statistical analyses is provided in the supplementary material. Angiogenic growth factor production G gvc 10 3 10 3 -10 5 Table 2. Parameter values of the factors describing the aberrant cellular behavior of NF1-mutated cells.