Pre-treatment growth and IGF-I deficiency as main predictors of response to growth hormone therapy in neural models.

Mathematical models have been applied in prediction of growth hormone treatment effectiveness in children since the end of 1990s. Usually they were multiple linear regression models; however, there are also examples derived by empirical non-linear methods. Proposed solution consists in application of machine learning technique - artificial neural networks - to analyse this problem. This new methodology, contrary to previous ones, allows detection of both linear and non-linear dependencies without assuming their character a priori The aims of this work included: development of models predicting separately growth during 1st year of treatment and final height as well as identification of important predictors and in-depth analysis of their influence on treatment's effectiveness. The models were derived on the basis of clinical data of 272 patients treated for at least 1 year, 133 of whom have already attained final height. Starting from models containing 17 and 20 potential predictors, respectively for 1st year and final height model, we were able to reduce their number to 9 and 10. Basing on the final models, IGF-I concentration and earlier growth were indicated as belonging to most important predictors of response to GH therapy, while results of GH secretion tests were automatically excluded as insignificant. Moreover, majority of the dependencies were observed to be non-linear, thus using neural networks seems to be reasonable approach despite it being more complex than previously applied methods.

Introduction

Growth hormone (GH) therapy is widely approved in children with short stature caused by GH deficiency (GHD). It has also been documented that GH therapy may be beneficial in children with idiopathic short stature (ISS); however, significant variability in its effectiveness in different patients is observed (1, 2). Taking into account the inconsistent and sometimes disappointing effects of treatment, the need to create prediction models of growth response to GH therapy has been defined by Ranke as a new paradigm for GH treatment in 21st century (3). Moreover, the use of growth prediction models has been proposed as a personalized approach to GH treatment in clinical practice (4). Development of such models should allow to determine more precisely the factors that influence responsiveness to GH, to facilitate realistic expectations concerning the therapy outcome and to identify poor responders, as well as to individualize and optimize treatment modalities (3, 4, 5). Prediction of GH therapy effectiveness should be a part of strategy for dealing with children with short stature, subjected to GH therapy (5).

In recent years, a number of prediction models have been created for different groups of patients and for different time horizons. Most frequently the models have been derived for prepubertal children with GHD (6, 7, 8, 9, 10, 11), while others were developed either for children in different pubertal stages (12, 13) or separately for prepubertal and pubertal ones (11). In some models, first-year response to treatment has been the only output variable (7, 13), while in others, the predictions have been done for few initial years of therapy (6, 8, 9, 10) or for the attained FH (11, 12). In one model, the total pubertal growth has been predicted (14).

Previously published models have been included as input variables – either only the data available before GH therapy onset or also the information gathered during treatment (6, 10, 11, 13, 14). Besides auxological parameters, GH peak after pharmacological stimulation has been considered as potential predictor of the therapy outcome in all the above-mentioned models (6, 7, 8, 9, 10, 11, 12, 13, 14), in some of them also concentrations of IGF-I (7, 8, 9, 13) and IGFBP-3 (7, 9) have been processed.

The most frequently used method of deriving models of GH therapy outcomes has been multiple linear regression (MLR) (6, 7, 10, 11, 12, 13, 14); however, in two models, non-linear technique of empirical curve fitting has been applied (8, 9). The main limitation of MLR models is the necessity to assume linear dependencies between input and output variables, while the real relationships between them can be non-linear and sometimes remain poorly defined or even unknown.

Development of advanced computational statistical methods (machine learning) seems to give the chance to overcome these limitations. Our research group has recently presented for the first time models of prediction of FH of GH-deficient children with the use of artificial neural networks (ANN) (15). Neural networks are very complex computational systems, designed to some extent to resemble neuronal connections in the brain. There are different types of ANN used for modeling of biological processes. In present study, multilayer perceptron (MLP) was selected as a model for prediction of response to GH treatment. This approach was chosen, because it does not require assuming any particular character of dependencies between explained variable (output) and explanatory ones (input). A neural model typically consists of:

input layer to which explanatory data are delivered,

hidden layers (1 or 2) that process the data,

output layer that delivers final result of prediction.

Each of these layers consists of several basic computational units – artificial neurons that process the data in order to produce desired result. In MLP hidden and output neurons calculate the weighted sum of all their inputs and transform it by sigmoid function. Every layer passes results of its calculations to the next until the output is reached and the solution presented to model’s user. The process of derivation of neural models, called training, consists of adjusting the weight coefficients of all neurons so that they produce the result that is as close as possible to desired one. Several algorithms exist that can be applied for this purpose. Before application, each model should be tested on a new data set, not used for its derivation.

From the practical point of view, the most important measures of prediction models accuracy are prediction errors (i.e. differences between real and predicted values of the output variable) and the proportion of variability explained and unexplained by the model. Very important is also the knowledge, which of the potential predictors have a real influence on the output value and which are redundant. Next essential issue is to determine the relationships between particular input variables and the output prediction. Finally, the model should enable clinicians to identify non-responders, including both the patients in whom the prediction is poor, as well as those, who do not improve growth despite good prediction.

Subjects and methods

Prediction models of GH therapy effectiveness were derived on the basis of data collected from 272 patients (197 boys, 75 girls) treated with GH due to isolated, non-acquired GHD in 20042016 for at least one year (1st year response model), including 133 patients (89 boys, 44 girls) treated up to the attainment of FH (final height model). In both models, most of children were prepubertal at therapy onset; however, 1st-year response model included 59 children who entered puberty before treatment (37 boys, 22 girls), while FH model included 39 ones (24 boys, 15 girls). The exclusion criteria were chronic diseases, multiple pituitary hormone deficiency, genetic syndromes, malnutrition, acquired GHD (brain tumors, cranial irradiation, injuries, etc.). Studies on GH therapy effectiveness were approved by the Committee of Ethics of Scientific Research in Polish Mother’s Memorial Hospital – Research Institute in Lodz. Written informed consent has been obtained from each patient and parents or legal guardians after full explanation of the purpose and nature of all procedures used.

Characteristics of the cohort of patients

The following data concerning the therapy onset or patients’ history were collected and pre-processed if necessary:

height at therapy onset (H ), expressed as height SDS for age and gender according to Polish reference data (H SDS ) (16);

gender (G), transformed to numerical values 0 for boys and 1 for girls;

chronological age at therapy onset (CA);

height at least 6 months before treatment (H ) and H SDS , respectively (16);

H SDS V – rate of changes of H SDS before treatment

body mass (M) and its SDS calculated with respect to height (M SDS) (16);

heights of mother (H ) and father (H ), expressed as H SDS and H SDS, respectively;

patient’s bone age (BA), assessed according to Greulich-Pyle standards (17), expressed as BA/CA ratio;

pubertal stage (PUB), transformed to numerical values: 0 – for prepubertal children, 1 – for pubertal ones;

gestational age (GA) and birth weight (BW), expressed as BW SDS for GA, according to (18);

GH peak after falling asleep (GH ) – an obligatory screening procedure in Poland, and in two stimulation tests: with clonidinie 0.15 mg/m2 orally (GH ) and with glucagon 30 μg/kg i.m., not exceeding 1.0 mg (GH ); GHD was defined as GH peak below 10.0 ng/mL in all tests; however, some patients with higher GH peak in one test were included to the study, as ones with either neurosecretory dysfunction or with normal GH peak only in nocturnal test (performed before its implementation as an obligatory screening procedure in Poland). This may constitute a bias of the study; however, we decided to include these patients as it was previously reported that assessment of nocturnal GH secretion may not identify all the patients with GHD (1, 19, 20); concentrations of GH were measured with IMMULITE 2000 system;

IGF-I and IGFBP-3 serum concentrations, expressed as IGF-I SDS and IGFBP-3 SDS for age and gender, respectively, assuming the log-normal distribution of IGF-I and IGFBP-3 concentrations, according to (21); concentrations of IGF-I and IGFBP-3 were measured with IMMULITE 2000 system;

initial GH dose in mg per week per kg of body mass (D); patients were treated with different initial doses of GH from the recommended range, chosen by their physicians, according to their clinical experience.

Additionally, the following data were obtained during or after 1st year of treatment and included as input variables in FH model:

patient’s height after 1 year of rhGH administration (H1), expressed as H1 SDS (16) and H SDS increase during 1st year of treatment (H SDS V1), the output variable in 1st year response model, calculated according to the formula:

IGF-I and IGFBP-3 after at least 3 months of treatment (IGF-I1 and IGFBP-31), expressed as IGF-I1 SDS and IGFBP-31 SDS, respectively (21). Finally, the differences between SDS of concentration before and during treatment were calculated and included in the model (ΔIGF SDS and ΔIGFBP-3 SDS).

The attainment of FH was confirmed by patients’ BA of 16 years for girls and 18 years for boys (22). The output variable in FH model was FH SDS (16).

Statistical characteristics of the whole group of patients are collected in Table 1 (data of subgroup who attained FH are presented in brackets). For the purpose of model derivation, whole datasets were randomly divided into 3 subsets:

Table 1

Statistical characteristics of patients’ cohort for 1st-year response model and for FH model (in brackets).

	Mean	s.d.	Minimum	Maximum
H SDS 0	−2.80 (−2.68)	0.68 (0.50)	−6.37 (−4.00)	−1.59 (−1.61)
H SDS V 0	−0.22 (−0.23)	0.23 (0.22)	−0.87 (−0.61)	0.58 (0.58)
CA (years)	10.9 (12.3)	2.8 (1.5)	4.2 (4.3)	14.6 (14.6)
M SDS	0.26 (0.65)	1.41 (1.26)	−2.8 (−2.08)	6.31 (3.61)
H M SDS	−1.00 (−0.94)	0.94 (0.85)	−4.43 (−3.03)	2.44 (1.93)
H F SDS	−0.92 (−1.06)	1.13 (1.12)	−5.09 (−5.09)	3.12 (3.12)
BA/CA	0.75 (0.80)	0.13 (0.10)	0.34 (0.46)	1.02 (1.02)
BW SDS	−0.96 (−0.91)	1.00 (0.96)	−5.38 (−3.8.)	1.37 (1.23)
GA (weeks)	39.3 (39.4)	1.3 (1.3)	35.0 (35.0)	43.0 (43.0)
GH N (ng/mL)	6.1 (6.3)	3.2 (3.7)	0.0 (0.0)	26.9 (26.9)
GH K (ng/mL)	7.1 (7.4)	4.0 (4.9)	0.4 (0.4)	40.0 (40.0)
GH G (ng/mL)	5.5 (6.0)	3.2 (3.7)	0.1 (0.1)	21.7 (21.7)
IGF-I SDS	−1.91 (−1.94)	1.1 (1.11)	−5.31 (−5.31)	0.78 (0.78)
IGFBP-3 SDS	−0.36 (−0.41)	0.95 (1.00)	−4.38 (−2.77)	3.16 (3.16)
D (mg/kg/week)	0.23 (0.23)	0.03 (0.03)	0.1 (0.13)	0.37 (0.37)
H SDS V 1	0.66 (0.57)	0.35 (0.29)	−0.09 (−0.09)	1.83 (1.45)
Δ IGF SDS	(2.54)	(0.98)	(0.26)	(6.37)
Δ IGFBP-3 SDS	(0.87)	(0.93)	(−3.87)	(3.17)
FH SDS	(−1.04)	(0.73)	(−3.19)	(0.60)

training set: 60% of data (163 cases for 1st year and 79 for FH prediction),

validation set: 20% of data (54 cases for 1st year and 27 for FH prediction),

testing set: 20% of data (55 cases for 1st year and 27 for FH prediction).

Neural models derivation and analysis

MLP networks for both 1st-year response and FH prediction were trained either by backpropagation (23) or conjugate gradient (24) algorithm. Parameters of networks (number of hidden layers, number of neurons) and training were selected by genetic algorithm (25), that consists of following sequence of operations:

Generation of initial population of networks (in our case 200 of them) with parameters set randomly from the range of allowed values.

Training and evaluation of obtained models by calculation of root mean square error (RMSE).

Elimination of certain portion of the least accurate networks (in this work 30%).

Generation of equal number of new models whose parameters are inherited from pairs of remaining networks (each parameter is assigned randomly to have a value coming from one of the two ‘parent’ networks) and their training.

Repeating steps 3 and 4 several times (here 5 times).

The best networks, in terms of accuracy of prediction, were chosen from final population predicting 1st-year response to GH treatment and FH SDS. Those networks were subjected to input data reduction by input cancellation algorithm (26). It is a model-based approach of redundant explanatory variables removal that is realized in the following way:

RMSE is calculated for original dataset.

Each input variable is separately replaced in whole dataset with its mean calculated over all cases from this set.

RMSE is calculated for all modified dataset.

A variable is removed that causes the lowest increase in RMSE when replaced by its mean.

Above operations are repeated until only one variable remains in the model.

The assumption behind the algorithm is that changing the value of important predictor influences strongly the prediction and aggravates its accuracy. At the same time, an insignificant variable is not used by the model in producing the result, therefore, assigning any value (e.g. average over whole set) in its place does not change prediction much. The last step of using this method is selection of number of significant inputs. It is done on the basis of sequence of RMSE values (or coefficient of determination, R 2), calculated after removal of each variable. A rapid decrease in R 2 (or increase in RMSE) means that an important predictor was eliminated and indicates the moment when reduction should be finished. The process of redundant variables elimination is performed automatically and the whole procedure is similar; however – from the technical point of view – much complicated than step-by-step backward elimination of insignificant variables in MLR models.

The networks derived by input cancellation were investigated further in order to identify the role of each input in prediction of response to GH treatment. The task was realized by sensitivity analysis in a similar way to approaches presented in the studies on ecological predictions (27, 28). It consisted in plotting the relationships between input and predicted output (P HSDS V for 1st-year response and P FH SDS for final height) for whole range of values of each continuous explanatory variable, when all other such inputs assumed the mean value over whole dataset. Binary inputs (gender and pubertal stage) were assigned one of the possible values that are indicated for each figure.

Results

From the population of 200 networks for 1st year prediction, produced in final iteration of genetic algorithm, a network with 2 hidden units – MLP 17:17-2-1:1 (number of input variables: #n of input neurons – #n of hidden units in consecutive layers – #n of output neurons: #n of outputs) was selected as the most accurate one. Similar procedure developing a model for final height indicated MLP 20:20-3-2-1:1 as the best model with all potential predictors. MLP 17:17-2-1:1 was trained by backpropagation, while MLP 20:20-3-2-1:1 by conjugate gradients algorithm. Both of those networks are presented schematically in Fig. 1; input variables abbreviations are explained in the'Characteristics of the cohort of patients' section where they are indicated by bold font.

Figure 1

Schematic representations of models for 1st-year response to GH treatment (A) and final height (B) with all available variables.

Those two models were subjected to input reduction by input cancellation algorithm. Results of this operation are presented in Fig. 2 as values of coefficient of determination (R 2) after removal of each input. Basing on the plot for 1st year (Fig. 2A), it was observed that 7 variables can be eliminated with only slight decrease in testing R 2; therefore, it was decided that final model should contain 9 inputs (MLP 9:9-2-1:1). In the model for final height, 10 explaining variables could be removed without significant change to R 2 (Fig. 2B), so final model was MLP 10:10-3-2-1:1. Variables that remained in both models are illustrated in Fig. 2C, where blue ellipse indicates 1st-year model and pink one FH model. Both reduced networks are presented in Fig. 3. Summary of errors and coefficients of determination calculated for both full and selected reduced models is presented in Table 2.

Figure 2

Results of input cancellation. Panels A and B present R 2 coefficient (y axis) for model after removal of consequent inputs (x axis) presents number of removed inputs. (A) Result for 1st-year model, (B) Result for final height model. Order of input elimination is indicated above curves. Panel C shows variables that were included in reduced models; blue ellipse indicates 1st-year model, pink one final height model, the intersection of ellipses contains variables that are common for both models.

Figure 3

Schematic representations of reduced models for 1st-year response to rhGH treatment (A) and final height (B).

Table 2

RMSE and R 2 for full and reduced models predicting 1st-year response to GH treatment and FH SDS.

Prediction horizon	Network architecture	RMSE			R 2 (%)
		TR	V	TS	TR	V	TS
1st year	MLP 17:17-2-1:1	0.258	0.255	0.267	42.8	47.3	48.7
	MLP 9:9-2-1:1	0.254	0.249	0.277	44.5	49.7	45.0
FH	MLP 20:20-3-2-1:1	0.489	0.498	0.493	53.7	59.3	44.0
	MLP 10:10-3-2-1:1	0.485	0.476	0.498	54.5	62.9	44.9

MLP number (n#) of inputs: n# of input neurons, n# of neurons in consecutive hidden layers, n# of output neurons: n# of outputs.

TR, training set; TS, testing set; V, validation set.

After input reduction, a sensitivity analysis was performed for reduced models MLP 9:9-2-1:1 and MLP 10:10-3-2-1:1. Input–output dependencies in 1st-year model are presented in Fig. 4; all of them were done for prepubertal children, separately for girls (red curves) and boys (blue curves). The observed relationships are non-linear and slightly different for boys and girls.

Figure 4

Sensitivity analysis for 1st-year model – MLP 9:9-2-1:1.

Analogical analysis for FH model (MLP 10:10-3-2-1:1) is shown in Fig. 5. In this case, gender and pubertal stage were not included in the model as redundant variables; therefore, the results are valid for both boys and girls either prepubertal or pubertal.

Figure 5

Sensitivity analysis for FH model – MLP 10:10-3-2-1:1.

Finally, we were interested in accuracy of prediction for particular patients. It was analyzed with the use of Bland–Altman plots (Fig. 6). For 95% of patients, prediction error did not exceed 0.5 s.d. in the model of 1st-year response to GH therapy and 1.0 s.d. in the model for FH; the mean error was about 0.27 s.d. and 0.50 s.d., respectively (for details see Table 2).

Figure 6

Bland-Altman plots for 1st-year model (A) and FH SDS model (B). Horizontal axis presents mean of real and predicted values, while vertical one prediction error (difference between real and predicted value). Green horizontal line indicates mean difference and red horizontal lines 95% confidence interval.

Discussion

In the present study, ANN model of 1st-year response to GH therapy explained 45% of its variability, with RMSE about 0.27 s.d./year that corresponds to about 1 cm/year. Other models of 1st-year response to treatment were created only for prepubertal children and explained from 33% (10) to 70% (7) of its variability, with the error of 1.46 cm (6) to 2.5 cm (7) in linear models and 0.23 s.d. in non-linear model (8). We did not manage to find any published model of 1st-year response to GH therapy for pubertal children or for children with different stages of puberty. Unfortunately, in some of the published models, the data on either explained variability or prediction error are missing. Moreover, despite the fact that the common model has been created for boys and girls, some gender-related differences in the relationships among the analyzed variables could be observed.

The best ANN model of FH SDS presented in the current study also explained 45% of its variability for testing group, with RMSE not exceeding 0.50 s.d. (that corresponds to about 3-cm difference between predicted and attained FH). For comparison, our previous model that included only pre-treatment data explained 42% variability of FH SDS with RMSE 0.63 s.d. Other models for FH prediction were derived either for only prepubertal children or separately for prepubertal and pubertal ones, which makes the direct comparison of our results with data from other studies impossible. However, the prediction error in our FH model was similar or even lower than that in published MLR models (11, 29). In our previous study (15), slightly better accuracy together with higher rate of explained variability of FH have been observed for ANN model vs MLR model, derived on the same data.

The best MLP model for 1st-year response contained 9 input variables selected out of 17 potential ones, as the most significant input variables there were selected: patients’ age at therapy onset, gender and pre-treatment IGF-I SDS, while the model eliminated all the input data concerning direct assessment of GH secretion. The best MLP model for FH prediction eliminated 10 out of 20 potential predictors. In this model, the most important variables were patients’ height at therapy onset (H SDS), pre-treatment IGF-I SDS and IGFBP-3 SDS, but again, none of GH stimulation tests results.

The problem of eliminating the variables that are considered as crucial for the diagnosis of GHD has been observed also in our previous ANN model for FH prediction, based on pre-treatment data only (15), as well as in some of the previously published linear and non-linear models (6, 11, 13, 14). Moreover, the same situation was reported in the model for patients with ISS (29) and in some models for patients with the broad spectrum of GH secretion capacities (8, 9). This finding requires further investigations on the significance of the results of GH stimulation tests for GH therapy effectiveness. It seems that the variability of GH response to stimulation should be taken into account, as previous studies have documented both early normalization of decreased GH secretion (30, 31, 32) and the converse situation (33).

In current classifications, GHD is defined as secondary IGF-I deficiency (34, 35). Nonetheless, IGF-I secretion was included as a potential predictor of GH therapy effectiveness only in a few models (7, 8, 9, 13, 14, 15). More severe IGF-I deficiency presented to be a significant predictor of better response to treatment in all the models except ones with incomplete data on IGF-I concentrations (8, 9). However, in the present study, IGF-I deficiency was not required as inclusion criterion, as in Poland, it is still not necessary to document decreased IGF-I secretion to diagnose GHD. Our results confirm the significance of IGF-I deficiency for better effectiveness of GH therapy. The importance of IGFBP-3 assessment in diagnosing GHD is still a matter of discussion (1, 36). This parameter presented to be a significant input variable in our FH model but not in 1st-year response one. In other studies, IGFBP-3 was included only in two models (7, 9), but it presented to be insignificant.

Creating FH model, we have entered as potential predictors the data obtained during 1st-year of treatment: an increase of H SDS (H SDS V ), IGF-I SDS (Δ IGF SDS) and IGFBP-3 SDS (Δ IGFBP-3 SDS); however, only H SDS V has been included in the final model. With respect to IGF-I significance, similar findings have been presented by Schonau and coworkers. (7). In all published models in which an increase of growth rate during GH therapy in the time period preceding the forecast was included, it has been the most important input variable (6, 10, 11, 29) or even the only significant predictor (9).

The complexity of neural models can seem to be an obstacle in their daily use in clinical practice. However, they could be implemented as user-friendly software as it was done for instance in the case of the models presented by Elias and coworkers (37). Still, application of models in the field of GH treatment response prediction requires some caution due to usage of various methods of hormonal assessment and differences in growth between populations. Certainly, any model should require calibration before it can be applied to children of different ethnicity or examined by laboratory methods with different reference ranges than those used in patients whose data we used to derive that model.

The undoubted contribution of our study to the deepening knowledge of the importance of particular auxological indices and results of hormonal tests for GH therapy effectiveness are presentations of non-linear dependencies between the analyzed input and output variables. This result is clearly a consequence of application of neural models that are more flexible in finding dependencies in data than other computational methods. Moreover, the progress of information technologies observed in recent years allows relatively easy translation of such complex tools to clinical practice.

Declaration of interest

The authors declare that there is no conflict of interest that could be perceived as prejudicing the impartiality of the research reported.

Funding

The study was partially supported by statutory funds from Polish Mother’s Memorial Hospital Research Institute, Lodz, Poland.