Proposed methodology for estimating w from measurements

In order to address the aforementioned ill-posedness, the following assumptions are made to reduce the dimension of the estimation problem. First, a L2-regularized least-squares problem is solved to estimate w by handling potential overfitting issues [37–39]. Second, this study aims to infer the values of w at the time instants only when the measurements are available (i.e., t_s). Third, instead of adding correction terms to all the model states, correction terms are given to only a subset of states, which will be denoted as $$ {\mathit{x}}_c\in {\mathbb{R}}^{n_s}$$ n_s < n_x, is selected a priori. Lastly, a linear interpolation is used to compute values of the correction terms at time t, when the measurements are not available (i.e., t ∉ t_s), as follows:

With these assumptions, the estimation of w(t) is reformulated as follows:

In Eq 5, the optimal value of α is unknown beforehand, so its value is optimized by five-fold cross-validation. Specifically, the available measurements $$ \mathcal{Y}=\widehat{y}$$ are divided into training and validation datasets with a 8:2 ratio in five different ways, and the regularized least-square problem (Eq 5) is solved with one particular value of α with respect to each of the five training datasets. Then, the optimal value of α is chosen by examining average model errors with respect to both the training and validation datasets.

Selection of x_c

A key step before solving Eq 5 is to identify x_c, whose trajectories are corrected by w_c. Specifically, two questions need to be addressed: first, what is the dimension of x_c (i.e., n_s), and second, when the value of n_s is known, which states in x should be selected to form x_c. In this study, we employ the idea of invertibility and a graph-theoretical approach to determine x_c.

Specifically, we choose x_c from x so that the resulted system is close to be invertible [40]. If a system is invertible, for a given value of x₀, unique values of y will correspond to unique values of inputs, so one could reconstruct the values of inputs from available output measurements [41, 42]. If w_c in the hybrid model (Eqs 5b and 5c) is viewed as an input to the system and the hybrid model is invertible, the values of w_c can be uniquely characterized from given measurements [41, 43]. Hence, it is our best interest to select the dimension of w_c as well as its placement so that the resulted hybrid model is invertible, which will attenuate the ill-posedness of the inverse problem (i.e., Eq 5). In this regard, Daoutidis and Kravaris [41] have shown that a dynamic system is invertible when the following matrix is nonsingular:

In this study, x_c is chosen so that C(x) of Eqs 5b and 5c is nonsingular, which will minimize the ill-posedness of the inverse problem (Eq 5). First, we let the size of x_c equal to the size of outputs (i.e., n_s = n_y) so that C(x) is square. Second, an optimal choice of x_c is identified systematically from x. In this regard, this study will assess the optimality of x_c through the following criteria that are based on the relative order:

The value of $$ {\sum }_{i=1}^{n_s}{r}_i$$ is minimum.

The value of $$ {\sum }_i^{n_s}{\sum }_{j\ne i}{r}_i/r_{ij}$$ is minimum.

Previous studies have demonstrated that the relative order measures ‘physical closeness’ between a correction term and an output [43, 45]: a smaller value of r_i represents a stronger connection between w_c and y_i. So, the first criterion renders w_c to have the maximum correlations with the outputs. On the other hand, the second criterion is to render each w_c,j will have the strongest correlation with only one output while having the weakest correlation with the remaining outputs [45, 46]. Consequently, by meeting the above two criteria, each correction term in w_c will have the strongest correlation with only one correction term, which will minimize the likelihood of the ill-posedness of Eq 5 [47].

In the perspective of the invertibility, selecting x_c based on the above two criteria, particularly the second one, maximizes the likelihood of C(x) to be nonsingular. To understand this point more clearly, the outputs and x_c are re-ordered so that the smallest element in each row is diagonally located in the relative order matrix, and C(x) is also rearranged correspondingly. By minimizing the second criterion, the possibility that values of all r_ij, ∀j ≠ i, are larger than r_i is maximized; therefore, the non-diagonal entries in the rearranged C(x) are likely to be zero. As a result, C(x) will be close to be a square diagonal matrix. Thus, achieving the aforementioned two criteria would guarantee the hybrid model (Eqs 5b and 5c) to be invertible as well as the reliability of the solution to Eq 5.

In order to select a combination based on the above criteria, the following steps are taken:

Enumerate all n_y permutations of x.

For each candidate, construct the corresponding relative order matrix, and compute its $$ {\sum }_{i=1}^{n_y}{r}_i$$.

Compute their ∑_i ∑_j≠i r_i/r_ij values.

Find the candidate that satisfies the condition.

For implementing the above procedure, a relative order matrix has to be constructed for each candidate to calculate the two criteria. However, performing Lie differentiation can be computationally expensive; hence, a graphical approach is implemented to evaluate relative orders of a system [43], which will be discussed in the next section.

A graphical approach to evaluate relative orders

A state-space model of a process (Eqs 5b and 5c) can be represented by a directed graph, which is defined by a set of vertices and a set of edges by the following rules [43, 48]:

States ($$ \mathit{x}\in {\mathbb{R}}^{n_x}$$), outputs ($$ \mathit{y}\in {\mathbb{R}}^{n_y}$$), and manipulated inputs ($$ {\mathit{w}}_c\in {\mathbb{R}}^{n_s}$$) are represented by a set of vertices in a graph.

If ∂f_i(x)/∂x_j ≠ 0, i, j = 1, …, n_x, there is a unidirectional edge pointing from the vertex of x_j to that of x_i.

If ∂f_i(x)/∂w_c,k ≠ 0, k = 1, …, n_s, there is a unidirectional edge pointing from the vertex of w_c,k to that of x_i.

If ∂y_l/∂x_j ≠ 0, l = 1, …, n_y, there is a unidirectional edge pointing from the vertex of x_j to that of y_l.

A path from one vertex to another is a sequence of edges without repeating vertices, and the path length is the number of edges included in one particular path.

Previously, Daoutidis and Kravaris [43] demonstrated that r_ij can be calculated by computing the shortest path length from an input w_c,j to an output y_i as follows:

Fig 1

Schematic illustration of the directional graph of a state-space model.

In summary, the procedures for selecting x_c, whose dynamics will be corrected by w_c, are as follows:

Set n_s, the size of x_c, to be equal to the dimension of outputs (i.e., n_y).

Enumerate all n_s permutations of x as candidates for x_c.

Construct a directional graph by adding correction terms to each x_c candidate enumerated in the previous step.

Construct the corresponding relative order matrix based on Eq 10.

Find a configuration that has the lowest $$ {\sum }_{i=1}^{n_y}{r}_i$$ and ∑_i ∑_j≠i r_i/r_ij values.

It should be noted that the identification of x_c by the relative order and graph theory can also guide the future model refinement. Specifically, since the identified x_c has the highest correlations with the output, further literature survey and experimentation on these states can be implemented to improve the differential equations for these states and thus increase the overall prediction accuracy of the first-principle model.

Model selection and training

The goal of the ANN training is to estimate hyperparameters of an ANN model, which include α, b, β, and c in Eqs 12 and 13 from available datasets. Here, the datasets include ANN input and output datasets, where the ANN inputs refer to t and x(t) from simulating the original first-principle model (Eq 1) while the ANN outputs refer to w estimated from solving Eq 5. All the ANN training sessions in this study are performed in the MATLAB Neural Network Toolbox with the Levenberg-Marquardt algorithm.

Since the structure of an ANN (i.e., the numbers of hidden layers and neurons in each hidden layer) is unknown beforehand, a number of ANN models with different structures are trained and compared to find the best one through evaluating their average corrected Akaike information criterion (AIC_c) [57, 58]. For a model with p number of parameters, AIC_c can be calculated as follows:

To find an optimal structure, datasets are randomly partitioned into the training, testing, and validation sets with a 70:15:15 ratio 100 times, and ANN models with different structures are trained 100 times to compute their average AIC_c [59, 60]. Then, the ANN structure with the minimum AIC_c is selected as the optimal one.

Once an ANN is developed, the final mathematical form of a hybrid model can be described as follows:

In silico measurements

The TNFα signaling pathway represented in Fig 3a can be described by the following mathematical model [61, 62]:

Table 1

Nominal parameter values of the correct TNFα signaling model shown in Fig 3a.

Parameter	Value	Parameter	Value
a1	0.6	b1	0.4
a2	0.2	b2	0.7
a3	0.2	b3	0.3
a4	0.5	b4	0.5
		b5	0.4

In this case study, the model (Eq 16) is considered as the true system, and it is used to generate in silico experimental measurements. Specifically, the NFκB activity (i.e., x₃) is measured every hour from 0 to 14 hours under three different input conditions (i.e., u = 0.5, 1, 2). It should be noted that the value of u is fixed while generating the in silico measurements. Also, experimental noise is introduced as follows:

Available first-principle model

On the other hand, we assume that the system is only partially understood, and Fig 3b represents the current understanding of the system, which is described by the following ODE model:

Fig 4

System-model mismatches in the TNFα signaling pathway model.

Comparison between experimental measurements (empty circles) and model predictions (solid lines) with 0.5 (blue), 1 (red), and 2 (black) of TNFα. The measurements are generated by simulating Eq 16 while the model predictions are obtained from Eq 19.

Hybrid model development

The first step of the proposed methodology is to determine a set of states whose trajectories need to be corrected by the addition of w_c. Since there is only one output (i.e., x₃), the number of states to be corrected by correction terms is one as described earlier. Then, a graphical approach is implemented to determine which state should be corrected by w_c. Fig 3c is a graphical representation of the first-principle model (Eq 19) to visualize the interconnections among the states. Since the number of the correction term to be added is determined to be one, the correction term can be added to one of x₁, x₂, and x₄ as shown in S1 Fig. In Fig 3, it is clear that there is only one directed edge pointing to the output (x₃), which stems from x₄. Consequently, the only feasible configuration for the correction term in this system is the first one in S1 Fig, where the correction term is placed to adjust the dynamics of x₄ and eventually the system output. It should be noted that the correction term is added to the differential equation of x₄ only.

Next, the regularized least-squares problem is solved to estimate the values of w_c at the time instants when the measurements are taken under three different input conditions. As the α value in Eq 5 for this system is unknown, its optimal value is found by the five-fold cross-validation. Table 2 shows the average normalized mean squared errors (MSE) between model predictions and the measurements for five different α values, and the optimal α value is determined to be one. Hence, the w_c estimation results corresponding to α = 1 are used for the subsequent analysis and ANN development. Before constructing an ANN, the accuracy of the estimated values of w_c is assessed by comparing the experimental measurements and the dynamics predicted by the available (incorrect) first-principle model (Eq 19) coupled with the estimated w_c. Fig 5 shows that the discrepancy between the predicted dynamics and the experimental measurements is diminished and thus validates the w_c estimation results.

Table 2

Comparison of the normalized MSE at different α values for the first case study.

α	normalized MSE
0.001	0.777
0.01	0.150
0.1	0.072
1	5.75 × 10−3
10	9.77 × 10−3

Fig 5

Validation of the accuracy of inferred w_c.

Comparison between predicted (red solid line) and measured (blue empty circle) NFκB dynamics when the activity of TNFα equals to (a) 0.5, (b) 1, and (c) 2. Blue dash lines represent the model predictions without the correction terms (w_c).

As the last step of the hybrid model construction, an ANN is developed to compute the w_c value. Here, inputs to the ANN are selected to be the values of states and input as well as the current time (i.e., [x₁ x₂ x₃ x₄ u t]), and the ANN output is the value of w_c.

The optimal number of hidden layers as well as the number of neurons in each hidden layer is to be optimized. As outlined earlier in the methodology, the AIC_c criterion is used to determine the ANN structure. To reduce the combinatorial complexity, the number of neurons in each hidden layer and the number of hidden layers are limited to ten and two, respectively. Then, each ANN is trained 100 times to compute the average AIC_c value with 100 different initial conditions for its hyperparameters, and the optimal ANN structure is determined by finding an ANN structure which results in the minimum average AIC_c value. Fig 6 plots the average AIC_c values for all possible ANN structures, and the AIC_c value reaches its minimum with eight and four neurons in the first and second hidden layers, respectively; hence, this particular structure is used for the subsequent analysis.

Fig 6

The average AIC_c values for different ANN structures for the first case study.

The filled circle represents the minimum average AIC_c value.

Among the 100 different ANNs with the optimal structure that have been trained in the previous step, the best ANN is chosen based on its R² statistic value. Specifically, an ANN with the highest R² value is chosen. As shown in Table 3, the R² statistics of the best ANN are sufficiently high to ensure its accuracy. Then, this ANN is coupled with the available (incorrect) first-principle model (Eq 19) to finalize the hybrid model development. In order to validate the prediction accuracy of this hybrid model, it is simulated under three input conditions and compared with the experimental measurements. As shown in Fig 7, the developed hybrid model can describe the true system dynamics fairly accurately under all the conditions, and the MSE of the model prediction is reduced to 0.00027 from 0.12 which is the MSE value of the original first-principle model (Eq 19). This result shows that the hybrid model constructed by the proposed methodology can accurately describe the system dynamics even when there is limited knowledge on the underlying system (i.e., only a model with model-system mismatches is available from the literature).

Table 3

The R² statistic values of the best ANN.

Training dataset	Validation dataset	Test dataset	Overall dataset
0.997	0.994	0.987	0.994

Fig 7

Validation of the developed hybrid model for the first case study.

Comparison between experimental measurements (empty circles) and hybrid model predictions (solid lines) with 0.5 (blue), 1 (red), and 2 (black) of TNFα.

Prediction capability of the hybrid model

The hybrid model is analyzed further in this subsection to assess whether the developed hybrid model possess the desired features of a hybrid model: intrepretability and generalized prediction capability.

First, the intrepretability of the developed hybrid model is tested by simulating and examining the temporal profiles of unmeasured states. Specifically, the dynamics of x₁, x₂, and x₄ predicted by the developed hybrid model are compared with those of the true system (Eq 16) to assess whether the hybrid model can be used in predicting unmeasured states. As shown in Fig 8, the predictions from the developed hybrid model agree well with the true system dynamics (Eq 16) and show significant improvement in the prediction accuracy compared with the available (incorrect) first principle model (Eq 19). Particularly, it is remarkable to note that the hybrid model can predict the dynamics of x₁ and x₂ quite well even though w_c is added only to x₄ for correcting its dynamics. Such improvement is possible since the hybrid model has incorporated the existing (known) interactions among the model states via the first-principle model (Eq 19).

Fig 8

Validation of the interpretability of the hybrid model.

The developed hybrid model is used to infer the dynamics of unmeasured states (i.e., (a) x₁, (b) x₂, and (c) x₄) under two input conditions: u = 0.5 (red lines) and u = 2 (black lines). The hybrid model predictions are compared with the true system dynamics (Eq 16) and the available (incorrect) first-principle model (Eq 19).

Additionally, it is of great interest to know whether the developed hybrid model has a generalized prediction capability. To this end, the hybrid model is used to predict the dynamics under three different input conditions (i.e., u = 0.7, 1.3, 1.7), and the model predictions are compared with the true system dynamics obtained from Eq 16. Here, these three particular input conditions are chosen since they lie within the input range used for training the ANN, but they are not identical to the inputs. As shown in Fig 9, the first-principle model as expected fails to capture the NFκB dynamics accurately. However, the hybrid model is capable of predicting the dynamics of x₄ fairly accurately. These results show that the hybrid model possesses the generalized prediction capability due to the incorporation of an ANN. In summary, this case study highlights advantages of using a hybrid model over a purely data-driven model or first-principle one since a hybrid model can essentially improve the overall model prediction accuracy while maintaining the model interpretability.

Fig 9

Generalized prediction capability of the developed hybrid model for the first case study.

The prediction accuracy of the hybrid model is assessed by comparing with true system dynamics when the activity of TNFα equals to (a) 0.7 and (b) 1.7. The dynamics predicted by the original first-principle model (Eq 19) are also plotted for the comparison.

Robustness of hybrid model

Additionally, since it is imperative to understand how the presence of measurement noise might impact the performance of the proposed methodology, different levels of noise are introduced by varying the distributions for μ_× and μ₊ in Eq 18 to examine the impact of the measurement noise. Specifically, a noiseless dataset (i.e., μ_× = 1 and μ₊ = 0) and that with a higher noise level ($$ ln{\mu }_{\times }\sim \mathcal{N}(0,0.05)$$ and $$ {\mu }_{+}\sim \mathcal{N}(0,0.05)$$) are generated, and their corresponding hybrid models are constructed. The additional datasets as well as the original dataset are plotted in S2 Fig to illustrate their difference in the noise level. With these additional datasets, the proposed methodology is implemented to develop the corresponding hybrid models, and S3 Fig illustrates the model accuracy of the hybrid models developed based on the noiseless and more noisy datasets, respectively. MSE of the developed hybrid models based on the noiseless and more noisy measurements is 0.00015 and 0.0048, respectively, while that of the hybrid model in Fig 7 is 0.00027. Overall, regardless of the noise level in the measurements, the developed hybrid models had significantly improved prediction capability as their predicted dynamics show reasonable agreements with the measurements. To further assess the impacts of the measurement noise, the interpretability of the hybrid models is assessed by using them to predict the dynamics of unobserved states as shown in S4 and S5 Figs. Overall, all the hybrid models were able to predict the dynamics of the unobserved states with reasonable accuracy.

However, it should be noted that both in the MSE and intrepretability analysis, the oscillations in the predicted dynamics become more noticeable with a higher level of noise. Since the oscillations are not present in the true dynamics, this comparison shows that the noise level may negatively influence the accuracy of a hybrid model as well as its interpretability. Also, even the hybrid model developed from the noiseless measurements produces the dynamics with nontrivial oscillations, which indicates that the ANN might have been overfitted. In order to mitigate such a problem, the future work in this direction could incorporate the following ideas to improve the hybrid model performance. First, a de-noise technique can be implemented prior to the hybrid model development so that the noise in the measurements will become less influential. In the literature, various methods such as finite difference with polynomial spline [67], spectral transformation [68], sparse Bayesian regression [69], and neural networks [70] have been proposed. Second, alternative ANN training mechanisms can be implemented to improve the performance of an ANN. So far, an ANN is trained with the Levenberg-Marquardt algorithm through Matlab Neural Network Toolbox. Since the presence of the oscillations suggests the potential overfitting issues in the ANN training, Bayesian regularization training method [71] may be implemented to mitigate this issue.

In a systems biology study, the parameter estimation is usually performed to quantitatively calibrate an uncertain model based on available measurements. In the literature, numerous studies have proposed a wide variety of methods to efficiently estimate parameter values, and these methods have been implemented to successfully model various biological systems [2, 72, 73]. However, if there is significant model-system mismatch, the parameter estimation alone may not be enough for calibrating a model. In such a circumstance, the proposed hybrid modeling approach can be implemented. As an example, the parameter estimation is performed for this system to see whether it is enough to train the model (Eq 19). Here, the parameter estimation is preceded by global sensitivity analysis to determine a subset of parameters that are identifiable from available measurements (see [74, 75] for details). The result of the sensitivity analysis is shown in Table 4, where ϕ_D represents the sensitivity index of a parameter set. Based on the magnitude of the sensitivity index, it is determined that b₁ and b₅ are identifiable. Then, a least-squares problem is solved to estimate their values by minimizing the difference between the model predictions and the measurements (i.e., in silico measurements simulated from Eq 16). The estimated values of b₁ and b₅ are 1 and 0.30, respectively, and the accuracy of the calibrated model is assessed by comparing the model predictions with the measurements as show in S6 Fig. It is clear that the parameter estimation alone is not enough to overcome the underlying structural mismatch between the model and the true system for this signaling pathway. In summary, the parameter estimation alone may not be enough for quantitatively calibrating the model if there is a significant degree of model-system mismatch. And, the proposed hybrid modeling approach will be a valuable option to construct an accurate and physically meaningful model.

Table 4

The sensitivity analysis results.

Parameter subsets	ϕD
b1	17.25
b5	15.18
b1, b5	137.4

w_c estimation

As outlined in the previous section, the first step in developing a hybrid model is to determine the number of correction terms needed as well as to which states these correction terms should be added. As the number of outputs for the system of interest is two (i.e., TNFα and IκBα), the dimension of x_c is also two.

Second, all possible permutations of choosing two from 49 model states are enumerated, and their corresponding digraphs are constructed to compute their corresponding relative order matrices. Based on the constructed relative order matrices, the minimum value of $$ {\sum }_{i=1}^{n_y}{r}_i$$ is found to be two, and Table 5 lists all the configurations whose $$ {\sum }_{i=1}^{n_y}{r}_i$$ values are equal to two. It should be note that r_1j and r_2j compute the relative order with respect to IκBα and TNFα measurements, respectively, in Table 5. Based on the result presented in Table 5, x₅ and x₃₇ are chosen as the best candidates to add w_c since this pair has the lowest ∑_i ∑_j≠i r_i/r_ij value. Here, x₅ and x₃₇ represent concentrations of IκBα and TNFα transcripts, respectively. Therefore, adding correction terms to these states is a reasonable choice since the predicted dynamics of a protein will become more accurate with more understanding of its transcript dynamics.

Table 5

Selection of optimal w_c locations by minimizing relative-order based criteria.

xs1	xs2	r11	r12	r21	r22	$$	∑i ∑j≠i ri/rij
5	37	1	6	6	1	2	0.33
1	37	1	6	5	1	2	0.37
3	37	1	6	5	1	2	0.37
34	37	1	6	5	1	2	0.37
2	37	1	6	4	1	2	0.417
4	37	1	6	4	1	2	0.417
5	28	1	4	6	1	2	0.417
1	28	1	4	5	1	2	0.45
3	28	1	4	5	1	2	0.45
34	28	1	4	5	1	2	0.45
2	28	1	4	4	1	2	0.50
4	28	1	4	4	1	2	0.5

With x_c = [x₁ x₃₇], the regularized least-squares problem (Eq 5) is solved to estimate W that contains the values of w_c at eight time points for each LPS concentration. Since the value of the regularization parameter α in Eq 5 is unknown beforehand, its optimal value is determined by the five-fold cross-validation. Table 6 shows the values of normalized MSE with respect to different values of α. From this result, the optimal value of α is determined to be 0.001, and the estimated values of W obtained with α = 0.001 are considered to develop a hybrid model.

Table 6

Comparison of the normalized MSE at different α values for the second case study.

α	normalized MSE
1 × 10−5	0.0275
0.0001	0.0283
0.001	0.0181
0.01	0.0187
0.1	0.0235
1	0.0310
10	0.0390

Before constructing an ANN model, the accuracy of the inferred W is assessed by comparing the experimental measurements and predictions from the imperfect model coupled with the inferred W. It should be noted that w_c is linearly interpolated to obtain its values at the time instants when the measurements are not taken. Fig 11 compares the predicted and measured TNFα and IκBα dynamics under three different LPS concentrations. Additionally, the predictions of the model coupled with the inferred W are compared with those of the imperfect model without W. Fig 11 shows that the addition of W significantly improves the model accuracy across all three LPS concentrations. Specifically, the addition of W renders the model prediction to match with trends observed in the experiments, which show the sustained low concentration of IκBα. At the same time, the addition of W helps the model predictions agree much better with the measured TNFα dynamics. Overall, these comparisons have demonstrated that the integration of W greatly improves the predictive capability of the available (imperfect) first-principle model.

Fig 11

The model prediction accuracy is improved by coupling the available (imperfect) first-principle model with the estimated w_c.

Red solid lines and blue empty circles represent simulated and measured TNFα dynamics and IκBα, respectively, in the presence of BFA under three different LPS concentrations. Blue dash lines represent the model predictions without the correction terms (w_c).

Additionally, in order to demonstrate the necessity of choosing an optimal x_c, a different x_c is chosen, and their corresponding values of W are estimated. Specifically, x_c = {2, 17} is selected, and their selection criterion values are $$ {\sum }_{i=1}^{n_s}=4$$ and $$ {\sum }_i^{n_s}{\sum }_{j\ne i}{r}_i/r_{ij}=1.25$$, which are much higher than those of the optimal choice of x_c. Then, the first-principle model is coupled with the estimated w_c, and its predictions are compared with the available measurements. Fig 12 illustrates that addition of w_c does not significantly improve the accuracy of the model prediction, which highlights the validity of determining the optimal set of x_c from Table 5 based on the relative-order criteria.

Fig 12

Importance of selecting an optimal x_c.

w_c are added to x_c = {2, 17} as explained in the test. Red solid lines and blue empty circles represent simulated and measured TNFα dynamics and IκBα, respectively, in the presence of BFA under three LPS concentrations. Blue dash lines represent the model predictions without the correction terms (w_c).

ANN development

With the inferred W, an ANN is developed to generalize the relationship between the first-principle model and w_c values. Specifically, inputs and output of the ANN are [x(t), t] and w_c(t), respectively.

Next, the ANN structure is optimized by minimizing the average AIC_c values. Similar to the previous case study, we will limit the numbers of hidden layers and neurons in each hidden layer to two and ten, respectively, and each ANN structure is trained 100 times. Fig 13 plots the average AIC_c value for all possible ANN structures, and the one with three and six neurons in the first and second hidden layers, respectively, is shown to be optimal since it provides the minimal average AIC_c value. Then, among 100 different trained ANNs with this particular structure, the best ANN is selected by its R² statistics. Table 7 shows the R² values of the chosen ANN, which are all above 0.98 and thus demonstrate its prediction accuracy.

Fig 13

The average AIC_c values for different ANN structures.

The filled circle represents the minimum average AIC_c value.

Table 7

The R² statistic valued of the best ANN for the second case study.

Training dataset	Validation dataset	Test dataset	Overall dataset
0.998	0.999	0.986	0.997

The developed ANN is then coupled with the available (imperfect) first principle-model to finalize the hybrid model. Fig 14 compare the predicted dynamics from the resulted hybrid model with the experimental measurements. The normalized MSE values of the hybrid model are 1.95 and 6.3 × 10⁻⁴ with respect to the TNFα and IκB measurements, respectively, while those of the original first-principle model are 12.5 and 7.0 × 10⁻³, respectively.

Fig 14

Validation of the developed hybrid model.

The IκBα dynamics predicted from the hybrid model with the developed ANN (red solid line) are compared with the measurements (blue empty circle) in the presence of BFA under the LPS concentrations of (a) 10ng/mL, (b) 50ng/mL, and (c) 250ng/mL. The TNFα dynamics predicted from the hybrid model with the developed ANN (red solid line) are compared with the measurements (blue empty circle) under the LPS concentrations of (d) 10ng/mL, (e) 50ng/mL, and (f) 250ng/mL. Blue dash lines represent the model predictions without the ANN.

This shows that the hybrid model with the developed ANN has generalized the prediction capability of the first-principle model coupled with the experimentally inferred w_c; as a result, the developed hybrid model now can be utilized to predict the dynamics of the NFκB signaling pathway in new conditions.

Although the developed hybrid model greatly improves the prediction accuracy, the model-system mismatch still remains. Specifically, under LPS = 10 ng/mL, the predicted dynamics of both TNFα and IκBα do not perfectly agree with the measurements. Specifically, the hybrid model predicts a monotonic increase in the IκBα dynamics and thus attenuation of TNFα synthesis, which indicates that the NFκB activity gradually decays during this time period. On the other hand, the measurements show that the TNFα synthesis slows down beyond 240 minutes while the IκBα level decreases after 240 minutes, which is not consistent with the model predictions. There is an additional mismatch in the IκBα dynamics under LPS = 50 ng/mL. Specifically, the hybrid model predicts an oscillatory behavior with two troughs at 60 and 240 minutes and a peak at 120 minutes while the experimental measurements indicate a monotonic increase from 60 to 240 minutes without an intermediate peak. Overall, such remaining model-system mismatches demonstrate that the BFA addition has more pronounced effects on the overall signaling dynamics after 200 minutes. This has been documented in our previous work [75], where the dynamics of IκBα in the presence of BFA deviate from those of IκBα without BFA after 200 minutes. Increasing the dimension of w_c may further improve the prediction accuracy due the increase in the degree of freedom. Alternatively, the first-principle model can be modified further by incorporating known mechanisms of the BFA-induced signaling pathways to improve the first-principle model before estimating the values of w_c. Specifically, it was suggested that the addition of BFA will elicit another signaling pathway called ER stress signaling pathway [75], which will suppress the translation of IκBα. In the future, this mechanism can be incorporated into the first-principle model to improve the accuracy of the hybrid model.