Comparison of Piecewise Structural Equation Modeling and Bayesian Network for De Novo Construction of a Quantitative Adverse Outcome Pathway Network

can be cross-linked to form large AOP networks (AOPN). The strength of an AOP is determined by its evidence support for causalities between different events, and the stressor and biological domains of applicability. A well-developed AOP allows better utilization of data generated from cost-efficient new approach methodologies (NAMs), such as in vitro high-throughput screening, high-content omics, and computational models, to inform regulatory decision making, thereby reducing the number and cost for conventional toxicity assessment using laboratory animals.

can be cross-linked to form large AOP networks (AOPN).The strength of an AOP is determined by its evidence support for causalities between different events, and the stressor and biological domains of applicability.A well-developed AOP allows better utilization of data generated from cost-efficient new approach methodologies (NAMs), such as in vitro high-throughput screening, high-content omics, and computational models, to inform regulatory decision making, thereby reducing the number and cost for conventional toxicity assessment using laboratory animals.

Introduction
The adverse outcome pathway (AOP) framework was introduced to better utilize the existing (eco)toxicological information for mechanistically based risk assessment of environmental stressors (Ankley et al., 2010).An AOP is a pathway that links a cascade of causally related toxicological events, from a molecular initiating event (MIE), via multiple key events (KE) occurring at increasing levels of biological organization, to an comings of conventional SEM such as multivariate normality and independence of the observations (Pugesek et al., 2003).These key characteristics make PSEM a promising approach for network model construction.To date, no study has employed the PSEM approach for qAOPN development.
Probabilistic methods such as BN modeling are being increasingly used in ecotoxicology and risk assessment (Moe et al., 2021a;Maertens et al., 2022).BN models can also be represented by DAGs, where causal/empirical relationships and the associated uncertainty/variability are quantified in conditional probability tables (CPTs).BN inherently incorporate uncertainty and can integrate a variety of information types, including expert elicitation.A BN approach has been proposed as an alternative for quantification of AOPNs (Moe et al., 2021b).This approach is based on a combination of regression modeling and BN modeling and is typically less data-demanding than the systems biology models.The quantified AOP-BN model can be run in several directions: 1) prognostic inference, run forward from the stressor node to predict the AO level; 2) diagnostic inference, run backward from the AO node; and 3) omnidirectionally, run from the intermediate MIEs and/or KEs.The BN-qAOPN reported by Moe et al. (2021b) was able to predict the state of AO with high accuracy when run from intermediate KE nodes, but the model performance decreased when the starting point was further removed from the target node.
The present study used a previously published AOPN and empirical dataset to construct novel qAOPN models.Three research questions were addressed: 1) Which linear AOP in the network contributed the most to the AO? 2) Can any of the upstream KE (including MIE) accurately predict the AO? 3) What are the advantages and limitations of PSEM or BN in qAOPN development?The model construction workflows and outcomes of the two approaches, PSEM and BN, were compared to suggest the best practice for de novo qAOPN development.
A quantitative AOP network (qAOPN) is an advanced form of AOP that captures the quantitative relationships of the KEs and biological complexities in a causal network, enabling the prediction of regulatory relevant AO using a mathematical/ computational model.A number of mathematical and computational approaches have been proposed to facilitate the development of qAOPN models (reviewed in Perkins et al., 2019;Spinu et al., 2020), with regression-based deterministic approaches (Foran et al., 2019;Perkins et al., 2019;Song et al., 2020) and Bayesian network (BN)-based probabilistic approaches (Zgheib et al., 2019;Moe et al., 2021b;Song et al., 2021;Spinu et al., 2022) being extensively explored for de novo model construction.Other approaches such as systems biology models have also been introduced but may require substantial data support and mechanistic understanding of the biological processes (Zgheib et al., 2019).
Among these approaches, structural equation modeling (SEM) has been increasingly used in various types of studies as a means to infer causal processes (Grace et al., 2012).Various statistical techniques and tools have been used in the process of specifying and evaluating the models based on SEM (Kaplan, 2009).In traditional SEM, the relationships among variables (i.e., the linear coefficients) are estimated simultaneously in a single variance-covariance matrix.Although the approach is well developed, SEM can be computationally intensive (depending on the sizes of the variance-covariance matrix).Moreover, additional assumptions such as independence and normality of errors are generally violated in studies or hard to test when the sample size is small (Yuan and Bentler, 1998).Piecewise structural equation modeling (PSEM) is a more flexible and potentially more powerful technique of SEM, which was proposed as an alternate approach to traditional variance-covariance-based SEM in the early 2000s (Lefcheck, 2016).The PSEM approach is based on directed acyclic graphs (DAG) and not restricted by the short- to generate inferences about the entire network.Unlike traditional SEM, where the χ2 test is commonly used to compare the observed and fitted covariance matrices and evaluate the overall goodness-of-fit, the goodness-of-fit of a PSEM model is tested using the directed separation of a DAG, which tests whether the path coefficients are significantly different from zero or whether we are justified in excluding them (Shipley, 2000;Shipley and Douma, 2021).
The coefficients of the qAOPN were estimated following these steps: 1) The values of the variables were scaled by normalizing to the mean values in the control group (unexposed group) prior to analysis; 2) a random sample of the observed data was generated using the bootstrapping method by a random resampling with replacement (Varian, 2005); 3) the coefficient of each edge in the qAOPN (Fig. 1) was estimated based on the PSEM algorithm.Simple linear regression models were used for the edges connecting the continuous variables, and multivariable logistic regression models were used for the edges connecting KE3, KE6, and KE11 with AO; 4) the predictive ability in terms of area under the receiver operating characteristic (ROC) curve of the four linear AOPs (Fig. 1) was evaluated using logistic regression models.Since a negative effect of a node on the AO 2 Materials and methods

AOP network (AOPN)
An AOPN linking reactive oxygen species (ROS) production to mortality was used as the prototype for qAOPN development in this study.This AOPN has been submitted to the AOP repository database AOP-Wiki as AOP:327-330 1 .The AOPN contains four linear AOPs with a common MIE (ROS production) and AO (mortality, Fig. 1).The supporting empirical dataset for model construction was derived from a previously published laboratory experiment, where the model aquatic crustacean Daphnia magna (14-day-old females) was exposed to a gradient (0, 0.0008, 0.05, 0.1, 0.2, 0.4 w/m 2 ) of a well-known oxidative stressor, ultraviolet B (UVB) radiation, for 2 and 7 days (Song et al., 2020).The AOPN (Fig. 1) serves as a conceptual model and qualitative structure for both model construction approaches (PSEM and BN).

Piecewise structural equation modeling
In PSEM, each set of relationships is estimated independently (or locally).This process decomposes the qAOPN into the corresponding simple or multiple linear regressions for each response, each of which are evaluated separately and then later combined  (KER), quantify the relationship between predictor and response variables by regression modeling; 2) apply the fitted regression models with associated uncertainty to simulate a high number of response values along the predictor gradient; 3) use the simulated values to parameterize the conditional probability tables, which are the quantitative links (KERs) in the BN model (see Fig. 2 for a summary of the workflow).A general issue for the analysis of the KERs is to identify the appropriate data combinations for response-response relationships.For each KE in the dataset used here, the response variable was measured in at least three replicates (i.e., three individual organisms).However, since the test material must be processed for a given response variable, an individual organism could not be used to measure more than one response variable.Therefore, there was no obvious way to directly link the replicate of one KE to a given replicate of another KE.could be offset by a positive effect of another node through the pathway, or transmit to the next node, both additive (logit(AO) = Σβ i KE i + e) and multiplicative (logit(AO) = Πβ i KE i + e) models were adopted.These steps (Fig. 2) were repeated 1000 times, and 1000 estimations were obtained for each coefficient and area under the ROC curve (AUC).Finally, distribution, mean, and 95% confidence interval of the estimates were obtained based on the 1000 estimations.

Bayesian network modeling
The BN approach for constructing a qAOPN mainly followed the procedure described by Moe et al. (2021b) using a combination of explorative statistical modeling, parametric regression, and simulation to quantify all components of the BN model.The general steps are as follows: 1) For each key event relationship scribed above), two simpler methods were applied to quantify the BN for comparison.For both alternatives, the CPT of each KER was quantified with an equation specified directly in the BN analytical software Netica version 6.042 (Norsys Software Corp., Vancouver, Canada).The equation was based on the coefficients estimated by the GLM (as above) but without the estimated standard deviation.Instead, the uncertainty of each KER was generated by the built-in sampling uncertainty function in Netica.The first alternative BN, labelled EQU-LIN, had all KERs estimated as linear models except for those with mortality as the response (logistic models).This set of models corresponds to the models used for SIM-LIN and for PSEM.For the second alternative BN, labelled EQU_SEL, a more comprehensive model selection was performed for each KER.The candidate regression models were linear, log-linear, quadratic, and logistic (sigmoid).For each KER, the optimal model was selected based on the percentage of deviance explained, supported by visual inspection of the scatter plots and fitted GAM curves.
For nodes with more than one parent node, such as KE-6 and the AO, the parent nodes were assigned equal weights in the BN model.A differentiated weighting of multiple parent nodes, and thereby the different pathways of the qAOPN, can be further refined later.
Sensitivity analysis was applied to characterize the influence of individual parent nodes on the target node (adverse outcome -mortality).The sensitivity is measured as mutual information between the target node and the parent node, which corresponds to the reduction in entropy of the target node (measured in bits) due to a finding (evidence) at the parent node.
For example, replicate 1 of KE-4 may be linked to replicate 1, 2, or 3 of KE-5.The solution used here was to generate paired predictor-response (x-y) values for the KERs using all possible combinations of the replicates for each x-y relationship.For example, KER no.7 (KE-4 → KE-5) would have 3×3 = 9 data points.This approach inflated the number of data points but avoided the need for a subjective or random selection of x-y pairs.
After preprocessing the data, the following procedure was used for each KER individually: 1) The functional form of the x-y relationship was explored and characterized using non-parametric regression (generalized additive models -GAM).2) An appropriate parametric regression model based on inspection of the fitted GAM curve.A ranking of the alternative models by the percent of deviance explained compared to the null model was used to guide the model selection.To compare the BN approach with the PSEM, which is based on linear models, we considered only linear models or linear models with transformation of the predictor variable (quadratic, logit, and logarithmic), corresponding to generalized linear models (GLM).3) The estimated coefficients and standard error of the fitted regression model were used to simulate a large number of data points (10,000) along the x-axis.4) The x and y variables were discretized into 5-6 intervals required for implementation of BN.Regression tree analysis was used for identifying breakpoints in the x-y relationships as a basis for selecting the appropriate number of intervals and the interval boundaries for each variable.5) The frequency distribution of the simulated data points across the x-y grid cell were used to derive the CPT of the x-y relationship.
In addition to the method for quantification of the BN based on simulations from linear models (labelled SIM-LIN as de- gression models were fitted further, and the AUCs of both models were larger than 0.9 (Fig. 5), indicating strong predictive ability of AOP4.
Although one or two coefficients in AOP1, AOP2, and AOP3 were not statistically significant, they showed excellent predictive ability regarding the AO and all the AUCs were larger than 0.9 (Tab.1).

BN-qAOPN
The curves estimated by GAM together with the estimated standard errors (Fig. 6) were used to characterize the shape of the KERs.Thresholds estimated by regression tree modeling were applied to identify the most important breakpoints in the KER shape.
The selected node intervals guided by the regression tree analysis are shown as gridlines in Figure 6.Most of the KE nodes were discretized into 5 intervals, which was considered sufficient for capturing the KERs.The stressor node (UVB) and the MIE node (ROS production) were given a higher reso-

PSEM-qAOPN
By using the PSEM approach, the regression coefficients were estimated for all adjacent KEs, with the KEs in AOP3 and 4 being highly significant (p < 0.001) compared to those in AOP1 and 2 (Fig. 3).As all the coefficients between any two adjacent KEs in AOP4 were statistically significant, the prediction performance was further evaluated for this specific qAOP.Multiplicative and additive multivariable logistic regression models were fitted, and the AUCs of both models were larger than 0.9, indicating strong predictive ability of the qAOP.
The distributions and average values of the coefficients of the edges between the adjacent KEs and AO along the pathways in the qAOPN are shown in Figure 3.
Both multiplicative and additive multivariable logistic re-  or probability distributions for the subsequent nodes.For each child node, the prior probability distribution was determined by the conditional probability table, represented by an arrow from its parent node(s).When the BN was instantiated for the highest UV level (Fig. 7B), the probability distribution of the MIE (ROS production) shifted towards higher intervals (average increased from 19 to 30).The probability distributions of the subsequent KEs correspondingly shifted towards higher or lower values; the signal was fading for each step.In consequence, the predicted AO (mortality) was practically unaltered (average from 0.406 to 0.41).However, the BN can also be instantiated from one of the KEs, which are fewer steps away from the AO and therefore more of the signal is retained.For example, a situation with the highest interval of KE-6 (ATP, representing intracellular energy transfer) predicts a lower probability of mortality (average 0.35).Running the BN model from an intermediate node can also provide a prediction from lution (8 and 6 intervals, respectively) to better capture changes at the start of the AOP.The AO node was discretized into 5 equidistant intervals for easier interpretation of the outcome.
For the BN presented here, we attempted to align the structure with the SEM as closely as possible to facilitate comparison.Therefore, the conditional probability tables (CPTs) of all KEs were derived from linear models (y = b0 + b1* x), while the CPTs of the AO nodes were derived from logistic models (y = exp(b0 + b1* x)/(1 + exp(b0 + b1* x)).The estimated parameters are reported in Table 2.
The resulting structure of the BN model with the initial settings is shown in Figure 7A-C.The circles represent intermediate versions of the subsequent nodes (KE-6 and AO) belonging to different pathways.In both cases, the intermediate nodes are combined into their respective child nodes by equal weighting.Figure 7A shows the uninstantiated (default) version with uniform probability distribution for the stressor node and pri-  of uncertainty for the CPTs.Likewise, the influence of the individual KERs on the target node also depends on the method of BN parameterization.
Considering the sum of percent mutual information for all the nodes combined (not shown), the BN based on CPTs from simulations had higher sensitivity (69%) than the BN based on linear models with built-in sampling uncertainty (EQU-LIN, 45%), but slightly lower than the BN based on selected best-fitting regression models (EQU-SEL, 73%).

Discussion
The present study employed two directed acyclic graph (DAG) types of model construction approaches, piecewise structural equation modeling (PSEM) and Bayesian network (BN) for de novo construction of a complex quantitative adverse outcome pathway network (qAOPN).Both approaches were able to quantify the response-response relationships and assist the identification of the most influencing AOP in a network.The advantages and limitations of these approaches are discussed below.

Piecewise structural equation modeling
Although PSEM has been successfully applied in some studies, such as determination of the links between the environmental conditions and diatom diazotroph associations (Stenegren et al., 2017) and evaluation of the role of quantity implicature in child language (Grinstead et al., 2022), it is still new in the field of qAOP development.Traditional SEM requires multivariate normality and independence of the variables, like many other statistical tests.The new extension of the PSEM relaxes these assumptions, so it is now possible to fit models to a variety of non-normal distributions and model non-independence through specification of fixed correlation structures or even random ef-the parent nodes (diagnostic inference).For example, the highest interval of ATP corresponds to a lower stressor level than the initial setting: the probability of UVB exposure being below 0.075 increased from the initial 37.5% to 68%.
Two alternative versions of the BN model were also constructed (EQU-LIN and EQU-SEL) to explore the performance of the BN with different assumptions underlying the quantification of the CPTs (Tab.2).For all alternative parameterization, the sensitivity of the target node (the final AO node) to the stressor and to the MIE node were close to zero.An analysis of the model performance in terms of accuracy (correct predictions) when the model is run from these nodes is therefore not meaningful.Instead, the sensitivities of the different model components were inspected more closely.
For each pathway, the sensitivity score listed for an AOP represents the percentage mutual information between the intermediate AO node of the pathway (shown as grey circles in Fig. 7) and the target node (the final AO node).Likewise, the sensitivity listed for individual KE nodes is the mutual information between the KE node and the target node.In general, the target node had very low sensitivity to changes in the stressor node for all of the alternative parameterizations (Tab.2).This reflects the low percentage of deviance explained for some of the KERs in the network even when the relationships are estimated with non-parametric regression by GAM (Fig. 6), which in general provided a better fit than any of the parametric models.The comparison of the three alternative parameterizations shows that the rank of influence of the AOPs differs among the three alternatives.For example, when the KERs are forced to be linear, the pathway AOP1 has the highest influence; but when the best-fitting regression model is selected, the pathway AOP2+3 has higher influence.This means that the relative influence of the different AOPs in the network depends on the method used for estimation of coefficients and quantification

Tab. 2: Sensitivity analysis for the Bayesian network model with the three alternative parameterizations of conditional probability tables
The prefixes SIM-and EQU-refer to the source of uncertainty: SIM = from simulations based on regression model estimates; EQU = from the built-in sampling uncertainty function in Netica.The suffixes -LIN and -SEL refer to the model selection: LIN = linear; SEL = selected as best fitting.
KEs acted directly and in concert to influence the AO.In addition, no bidirectional relationships caused by a shared underlying driver have been considered, only crude approximate correlated errors are implemented, and cyclic relationships cannot be evaluated (Lefcheck, 2016).

Bayesian network
An inherent strength of BN models is the explicit and transparent handling of uncertainty in each step (node) of the model.In this study, for each KER, the predictor-response relationships as well as the variability was first explored by GAM and then quantified by generalized linear models (including linear regression).The uncertainty estimated by the selected GLM could then be incorporated into each KER of the BN via the CPT.In the AOP context, the BN represents a causal model, in a similar way as the PSEM.An additional benefit of the BN compared to the PSEM is the flexibility of the CPTs (causal links) to allow for more flexible relationships.Although the CPTs in this study were limited to GLMs for easier comparison with PSEM, a CPT can in principle represent any functional form.Moreover, a CPT can easily capture interactions (synergisms or antagonisms) between multiple parent nodes, although this aspect was not explored in the present study.fects (Lefcheck, 2016).Unlike traditional modeling approaches, PSEM is often preferred in a system where several predictors and outcomes are connected.It incorporates hypothesized causalities between the nodes in the network, which is based on previous knowledge of the investigated system, and it facilitates identification of the direct, indirect, and cascading effects in the system.Although PSEM is built on assumed causalities, whether the causalities are correct is difficult to test when the predictors and outcomes were observed at the same time.Fortunately, in the current experimental study, the predictors and outcome occurred sequentially, indicating that the causalities were not inversed.It should be noted that the PSEM also needs to meet all the assumptions of the underlying generalized linear models, such as homogeneity of variance, which was obtained in this study by scaling the variables during preprocessing of raw data.
While the PSEM is a considerable leap forward from traditional SEM, there are also limitations in the implementation of PSEM in the current study.First, only direct effects between the nodes were investigated due to the complexity of the network.Second, we did not incorporate latent or composite variables in the current study, which is indeed one of the great advantages of conventional SEM.Nevertheless, the qAOPs constructed by PSEM showed excellent predictive ability, suggesting that the  also of importance, both for the overall sensitivity of the BN and for the relative sensitivity of the individual AOPs in the network.An optimal workflow for qAOPN development would combine the more flexible modeling of KERs allowed by the BN with the higher predictive ability displayed by the PSEM.Apart from these, additional experimentation that captures the dose-response patterns across a time scale (i.e., temporal response) may be beneficial for constructing more robust causal models using these approaches.Although not achieved in the present study, validation of the qAOPN models by follow-up experimental studies would be highly appreciated in the future.

Conclusion
The present study has demonstrated that both piecewise structural equation modeling and Bayesian network are suitable for de novo construction of a complex qAOP network.Besides quantification of response-response relationships, these two approaches were also capable of identifying the most influencing linear AOP in a complex AOPN and evaluating the predictive ability of the qAOP.Future studies will also identify the most critical (influencing) KE(s) in the qAOPN and how accurately this KE can be used to predict the final adverse outcome that is relevant for hazard and risk assessment.The workflow described herein can be easily repeated by other laboratories for efficient qAOP model construction utilizing routine ecotoxicity data.
Another unique feature of BN models is the ability to run the model in multiple directions: prognostic (forward) from the stressor node or from the MIE; diagnostic (backward) from the AO, or omnidirectional from one or more intermediate KEs.This feature allows for more types of inferences from a qAOP(N) if the model has sufficiently high predictive ability.
The BN model developed here, however, did not obtain sufficient sensitivity to have predictive ability.There are several explanations for the lack of sensitivity.One reason is the low correlation of x-y values for some of the KERs (Fig. 6), resulting in a weak signal in the corresponding CPTs.A different handling of the replicates could have resulted in lower uncertainty in the BN.For example, randomly selecting only one replicate for each x and each y would potentially lead to lower variability and thereby better model fit.Alternatives for handling of response-response relationships were beyond the scope of this paper but will be addressed in later investigations.
Another obvious reason for low sensitivity of the BN model is the discretization of continuous variables, which inevitably results in loss of information and precision.The discretization of variables is a standard process and inherent shortcoming of traditional BN modeling.However, more advanced BN methods allow for continuous nodes (e.g., Jackson-Blake et al., 2022) or a combination of discrete and continuous nodes (e.g., Moe et al., 2020).The propagation of quantified uncertainty throughout the model is also a unique property of BN models, which provides transparency and reflects the modeler's best knowledge of the system but also contributes to low sensitivity of the target node (adverse outcome) to the earlier components of the network.
The comparison of the three methods for parameterization of CPTs (Tab.2) showed that the method with estimation of linear models with uncertainty and subsequent simulations (SIM-LIN) had overall higher sensitivity than the alternative with built-in uncertainty (EQU-LIN).This comparison suggests that the combined estimation and simulation method (SIM-LIN), which provides a precise characterization of the uncertainty in the linear KERs, has the potential to generate qAOPNs with high sensitivity, even though the resulting qAOPN in this example had low predictivity.However, the alternative method, where the best-fitting regression models were selected for generating CPTs (EQU-SEL), resulted in even higher sensitivity, suggesting that selecting appropriate non-linear parametric models for the KERs can be important for overall model sensitivity of the qAOPN.A repetition of this exercise on a new dataset with more observation and/or with less variation among replicates might result in a BN-AOP that also has higher predictivity.

Recommendations for future direction
Comparing the two approaches used in this paper, PSEM and BN, the PSEM obtained higher predictive ability than the BN for the given dataset.However, in spite of the low overall sensitivity of the BN model, some important properties were identified for the purpose of qAOPN development.The exploration of alternative methods for parameterization of the CPTs (Tab.3) showed that the selection of the parametric model underlying each CPT is

Fig. 2 :
Fig. 2: The workflows for de novo quantitative adverse outcome pathway network (qAOPN) construction using two parallel approaches, piecewise structural equation modeling (PSEM) and Bayesian network (BN) modeling CPT, conditional probability table; KER, key event relationship

Fig
Fig. 3: Distribution of the coefficients between the adjacent key events (KEs) and adverse outcome (AO) in the quantitative adverse outcome pathway network (qAOPN)

Fig. 4 :
Fig. 4: Average coefficients over the 1000 bootstrap iterations for the key event relationships

Fig. 5 :
Fig. 5: One of the 1000 evaluations of predictive ability of the quantitative adverse outcome pathway #4 (AOP4) in the network using multiplicative (A) and additive multivariable (B) logistic regression models

Fig. 6 :
Fig. 6: Observations (red dots) and fitted generalized additive models (GAM, black curves) with standard errors (SE, dashed blue lines) The numbers in the upper right corner are the percentage of deviance explained by the GAM model.

Fig. 7 :
Fig. 7: A Bayesian network (BN) representation of a quantitative adverse outcome pathway network (qAOPN) (A) The BN with prior probabilities and uniform probability distribution for the stressor node, UV. (B) The BN instantiated from the stressor node.(C) The BN instantiated from key event (KE) no.6, ATP.

Tab. 3 :
Parameter estimates from regression models from the piecewise structural equation modeling (PSEM) and Bayesian network (BN) workflowsNote that the PSEM used variables in normalized scale while the BN workflow used variables in original scale (raw data), which resulted in coefficients of different magnitudes for some of the key event relationships (KERs ). AO, adverse outcome; KE, key event; MIE, molecular initiating event