Integration of temporal single cell cellular stress response activity with logic-ODE modeling reveals activation of ATF4-CHOP axis as a critical predictor of drug-induced liver injury

Drug-induced liver injury (DILI) is the most prevalent adversity encountered in drug development and clinical settings leading to urgent needs to understand the underlying mechanisms. In this study, we have systematically investigated the dynamics of the activation of cellular stress response pathways and cell death outcomes upon exposure of a panel of liver toxicants using live cell imaging of fluorescent reporter cell lines. We established a comprehensive temporal dynamic response profile of a large set of BAC-GFP HepG2 cell lines representing the following components of stress signaling: i) unfolded protein response (UPR) [ATF4, XBP1, BIP and CHOP]; ii) oxidative stress [NRF2, SRXN1, HMOX1]; iii) DNA damage [P53, P21, BTG2, MDM2]; and iv) NF-κB pathway [A20, ICAM1]. We quantified the single cell GFP expression as a surrogate for endogenous protein expression using live cell imaging over > 60 h upon exposure to 14 DILI compounds at multiple concentrations. Using logicbased ordinary differential equation (Logic-ODE), we modelled the dynamic profiles of the different stress responses and extracted specific descriptors potentially predicting the progressive outcomes. We identified the activation of ATF4-CHOP axis of the UPR as the key pathway showing the highest correlation with cell death upon DILI compound perturbation. Knocking down main components of the UPR provided partial protection from compound-induced cytotoxicity, indicating a complex interplay among UPR components as well as other stress pathways. Our results suggest that a systematic analysis of the temporal dynamics of ATF4-CHOP axis activation can support the identification of DILI risk for new candidate drugs.


Introduction
Drug-induced liver injury (DILI) is one of the most frequently encountered drug-related toxicities in drug development and clinical usage [2,6]. Despite the higher incidence compared to other drug adverse reactions, the underlying molecular mechanisms by various drugs causing liver injury remain unclear due to multiple limitations such as complexity of the injuries, non-robust findings related to host-specific factors from the current in vitro approaches, and lack of correlation on liver injury markers between in vivo findings and clinically significant DILI [29]. Various key events are involved in DILI, including the release of toxic metabolites that damages neighboring non-parenchymal cells, activation of immune cells, formation of lipid droplets and release of chemokines to induce macrophage crowns [46]. On top of that, A major critical key event during liver toxicity is the onset of cell death of hepatocytes. The occurrence of cytotoxic events in DILI has been linked to the activation of specific cellular stress response pathways [34] including: inflammatory response (IR) pathways mediated by the NF-κB transcription factor [20,37], oxidative stress response (OSR) pathways mediated by nuclear translocation of NRF2 (NFE2L2) [9], DNA damage response (DDR) pathways mediated by TP53 activation [8,18], and the unfolded protein response (UPR) pathway mediated by ATF4, ATF6, and XBP1 transcriptional programs [10,21]. Understanding the dynamics and amplitude of the activation and interaction of these cellular stress response pathways in hepatic cells in relation to protection against or onset of cytotoxicity could provide better insights in critical determinants underlying DILI.
Numerous advanced in vitro test systems based on 2D and 3D (co)cultures have been employed to study DILI mechanisms [1,29]. Typically, experimental approaches with these test systems rely on endpoint measurement and do not capture the dynamics of cellular stress responses thereby limiting temporal mechanistic insights [14]. Our previous work has demonstrated the applicability of bacterial artificial chromosome (BAC)-based GFP HepG2 reporter cell lines to study the dynamics of cellular stress response activation [48][49][50]. This technology allows the expression of GFP-labeled proteins representing critical components of the different stress response pathways under the control of endogenous promoter activity, thus allowing high content imaging and quantification of the temporal protein expression at the single cell level using multiparametric automated image analysis pipelines [49]. These cellular stress response reporters have proven successful in gaining comprehensive insight in temporal dynamics of several specific stress response pathway components in the first twenty four hours after xenobiotic or cytokine exposure [16,37,50]. Here we applied the panel of reporters to uncover particular cellular stress response pathways that are of critical significance to predict liabilities for DILI-related cytotoxicity.
Concurrent with the development of in vitro models, recent mathematical modelling strategies have been developed to study the dynamics of signaling networks. The models capture details and offer mechanistic insights into molecular perturbations that link external signals, e.g. cellular stress to cellular outcomes [26]. Different areas of application of mathematical modelling include neurodegenerative diseases, cardiovascular diseases and cancer, which were demonstrated to have etiology based on the deregulation of signaling networks [19,25,31]. The investigation of deregulated signaling pathways in the field of pharmacology and toxicology with mathematical model is also gaining popularity and provides a better understanding on how signaling molecules are interacted among multiple pathways and how they are connected to liver-specific pathology which are observed at the cellular level [23,37,54]. Among the mathematical modelling approaches applied to model and analyze signaling networks, logic models offer a balance between required model inputs and computational time versus mechanistic insights in terms of network connectivity and their relationships which are depicted by logical gates AND, OR and NOT [32,52]. The simplest form of the logic model is a Boolean network, requiring a functional protein interaction network and discretized experimental data at quasi-steady-state as input in order to perform network inference and contextualization. Yet, this simple logic formulation is not suitable for time-course experimental data as generated by live cell imaging. To overcome such limitation, variants of logic models have been developed. One of the most advanced forms, the logic-based ordinary differential equations (Logic-ODE), derive a set of ODEs from logical rules and are capable of fitting quantitative time-course experimental data [51]. Compared to mechanistic-based ODE models which are derived from biochemical reactions, Logic-ODE models require no prior information on kinetic parameters and can be optimized with minimal parameterization. The optimization time and model calibration effort are also less complex for Logic-ODE models which make them suitable for large-scale modeling of quantitative time-course data with limited prior kinetic and biochemical information.

Table 1
List of the tested compounds. Table of the compounds tested in the experiment including the c-max, compound code used in this article, actual concentrations annotated from dose levels 1 to 8, and the compound classifications. Cmax was adapted from the previous publication [35,50]. Compound  Given the advantages of the Logic-ODE modelling framework and the capacity of BAC-GFP reporter cell lines to comprehensively depict the temporal activation dynamics of stress response pathways, we applied a combinatorial approach to delineate the regulatory mechanisms orchestrating the stress response signaling pathways activated during DILI episodes in relation to cytotoxicity. Quantification over time of the GFP activity of critical cellular stress signaling molecules combined with cell viability collected from time-lapses of our large panel of BAC-GFP HepG2 reporters encompassing 4 mains stress responses (oxidative stress, unfolded protein response, DNA damage response, and Inflammation response) were used to contextualize the model. The parameterization of our Logic-ODE model unraveled two aspects of DILI: i) regulation of stress response activation is compound dependent, and ii) the ATF4-CHOP axis of the UPR is the pathway that correlates most with DILI-related cytotoxicity. Small interference RNA-mediated perturbation of main components of the UPR pathway supported our mathematical model and indicated the central UPR transcription factor ATF6 as a counter-regulator of drug-induced cell death caused by DILI compounds.

Chemical and reagents
All chemicals were purchased from Sigma-Aldrich -The Netherlands; except for cisplatin (Ebewe -The Netherlands) and nefazodone (Sequoia Research Products -Pangbourne, United Kingdom). All compounds were dissolved in DMSO; except for mitomycin-C (DMSO-PBS), N-acetylcysteine (PBS) and for cisplatin which was already manufactured as a solution. TNFα was purchased from R&D System (Abingdon, United Kingdom). All compounds in DMSO were maintained as 500-fold stock such that the final exposure did not exceed 0.2% v/v DMSO. PowerUp SYBR green real time PCR master mix was purchased from ThermoFisher. All primers were purchased from Sigma-Aldrich -The Netherlands.

Cell lines
Human hepatoma (HepG2) cells were purchased from ATCC -Germany (clone HB8065) and maintained in DMEM high glucose (Fisher Scientific -Bleiswijk, The Netherland) supplemented with 10% (v/v) FBS (Fisher Scientific-Bleiswijk, The Netherlands), 250 U/ml penicillin and 25 µg/ml streptomycin (Fisher Scientific -Bleiswijk, The Netherlands) in humidified atmosphere at 37 degrees Celsius and 5% CO2 /air mixture. All the BAC-GFP HepG2 reporter cell lines were previously established and characterized [48]. The cells were used between passage 14 and 20. For live cell imaging, the cells were seeded in Greiner black µ-clear 384 well plates, at 8000 cells per well.

Compound exposure and live imaging
Two days after seeding in either 384 well plate or 96 well plate (for RNA interference experiments), the cells were stained with 100 ng/ml live Hoechst 33342 in complete DMEM high glucose overnight. At the day of exposure (three days after seeding), the medium containing Hoechst was refreshed with complete DMEM containing 0.2 µM propidium iodide (PI) and Annexin-V-Alexa 633 (AnV). The compounds were added as such that the final concentration was as it is shown in Table 1. The plates were thereafter imaged no longer than 1 h after exposure. For inflammatory reporters (HepG2 ICAM1 and HepG2 A20), TNFα final concentration at 10 ng/ml was added 8 h after the compound exposures and the plates were imaged no longer than 1 h after the TNFα addition. The plates were imaged for 65-72 h using a Nikon TiE2000 confocal laser microscope (laser: 647 nm, 540 nm, 488 nm, and 408 nm), equipped with automated stage and perfect focus system. During the imaging, the plates were maintained in humidified atmosphere at 37 degrees Celsius and 5% CO 2 /air mixture. The imaging was done with 20x magnification objective and performed every 1.5 h. Each plate only contained one reporter cell line exposed to the set of compounds and with 3 replicates per concentration.

Real time -quantitative PCR
RNA from exposed cells was isolated using the RNeasy kit (Qiagen) according to the manufacturer's protocol. cDNA synthesis was performed using the RevertAid H Minus First Strand cDNA synthesis kit from ThermoFisher Scientific. The RT-qPCR was performed with the mixture of cDNA and PowerUp SYBR green master mix (Thermo-Fisher) for 40 cycles. The RT-qPCR and the analysis was done with High Resolution Melt Module for QuantStudio 6 Flex Real-Time PCR system (Applied Biosystems) for quantifying the expression of DDIT3, HSPA5, XBP1, ATF4, ATF6, EIF2AK3, and GAPDH. GAPDH expression was used as the housekeeping gene control and the fold change is calculated with the 2 − ΔΔCt method. The primer sequences for the RT-qPCR are listed as followed:

Image analysis
The images were manually sorted to exclude images which do not fulfill the criteria for analysis: non-biological background (fluorescent fibers), loss of nuclear signal, and in-focus images. The quantitative image analysis was performed with ImageJ version 1.52p and CellProfiler version 2.2.0 [24]. Firstly, the nuclei per image were segmented watershed masked algorithm on ImageJ and thereafter processed with an in house developed CellProfiler module [48,53]. The results were stored as HDF5 files. Data analysis, quality control, and graphics were performed using the in-house developed R package h5CellProfiler. For each reporter, the nuclear Hoechst33342 intensity levels, GFP intensity, PI area, and Annexin V area were measured at the single cell level.

Data analysis
The GFP intensity from cell population means of each image was calculated based on the single cell results. In addition, for each plate, the GFP intensity of DMSO control (or DMEM for cisplatin), was calculated to determine the background. The fraction GFP positive cells was calculated per plate based on the cells which have GFP intensity higher than 2 times of mean value. For ICAM1 and A20, GFP positive cells were calculated based on the subtraction of the cells which showed GFP intensity higher than the thresholds with the cells showing GFP intensity lower than the threshold. For this purpose, the thresholds were derived from the GFP intensity value in the TNFα treatment ± 1 times of mean GFP intensity in the DMSO control. These values were thereafter normalized by subtracting them with the GFP positive cells from the TNFα treatment. This calculation allowed us to capture the decrease of the A20 and ICAM1 responses due to compound exposures. To quantify the fraction PI and AnV positive, the PI images and AnnexinV images were masked with 2 pixel dilated nuclei based on nuclear segmentation) to exclude the background staining noise. The area of each nucleus and the corresponding PI/AnV object were divided to obtain a PI and AnV/ nuclei ratio. PI/AnV positive for every cell were defined as a cell with >10% of PI/AnV area.
Normalization methods were applied to GFP intensity values and cell count. For the GFP intensity, the value was scaled between 0 and 1 with the formula: (X-Xmin_replicate)/(Xmax_replicate -Xmin_replicate, where X = GFP intensity value). For ICAM1 and A20, the value was scaled between − 1 to 1 to account for the up and down regulation of the TNFα-modulated ICAM1 and A20 regulation upon the exposure of the compounds. For cell count, the value was scaled from − 1 to 1 by calculating the log10 fold change of the cell count values from each observation to the cell count value of the control treatment (DMSO/ DMEM) at time point 0.
All features (normalized GFP intensity, fraction GFP positive cells, fraction PI positive cells, fraction AnV positive cells, and normalized cell count) were fitted by applying b-splines with 8 degree of freedom. After that, 60 discrete equidistant time points were selected to calculate means of replicates per time-point A hierarchical clustering was performed to define the compound cluster based on the temporal activation of stress response pathways and progressive cellular outcome (Annexin V signal, PI signal, and nuclear count) in each compound. First, we aggregated the value of every readout of each compound into one value by calculating the mean from the full concentration levels. Thereafter, euclidean-based distances between all time-course vectors were calculated. The mean euclidean-based distance from the all features was used as inputs for the WardD2-based clustering. AUC values were calculated from the time course results for each replicate, compound, and concentrations after the b-spline fitting with an in house Rscript. The correlation analysis was performed with Pearson correlation method with the chart.correlation function from the performance analytic package [38].

Statistical analysis
For all reporters and concentrations, three to four independent biological replicates from imaging experiments were performed. A twosample two-sided Student's t-test was applied to the AUC values from the time course results of the knocked-down experiment. Two-sided test was chosen accounting the responses upon the siRNA knock-down that can be either higher or lower than the control (siControl). The Student's t-test was performed to calculate the significance based on the (uncorrected) p-value between the siRNA target and siRNA control (siControl). The adjusted p-values thereafter were calculated with false discovery rate method.

Data representation
All high content imaging results were presented or modified with Illustrator CS6, ImageJ version 1.52p, and R (ggplot2 [47] and pheatmap [27]).

Mathematical modelling -data processing and normalization
The mean GFP signals derived from the integrated intensity channel within the cellular compartments which the respective molecules are active were background subtracted by the GFP signals from their solvents. All negative values which resulted from lower GFP signals from compound treatment compared to control were set to zero. The background subtracted GFP signals were then normalized to the maximal GFP signal per GFP reporter. Concentration of treated compounds were log10 transformed and scaled to the maximal values. The measurement for Apoptosis derived from Annexin-V signals was already measured in percentages and was directly applied as an input for modeling. All processed data are in a continuous range between 0 and 1 which is compatible with logic modeling.

Mathematical modelling -Network topology
The network topology of mathematical model was derived from the canonical pathways of the four stress response pathways connecting to all 13 measured GFP reporters [5,20,41,42]. Molecules in the canonical pathways which were not measured such as KEAP1 and ATF6 were excluded. 11 crosstalk interactions between the four stress response pathways were added from Omnipath [45], a collection of curated The signal intensities represent the responses from multiple proteins involved in particular stress response pathways. The cells were exposed to positive controls for each specific stress response (oxidative stress: DEM; unfolded protein response: tunicamycin; DNA damage response: etoposide; inflammatory response: diclofenac). Blue indicate nuclear staining (Hoechst) and green indicate GFP signal. The nuclear images for NRF2 and ATF4 are removed to enhance the clarity of GFP responses. B) The temporal dynamic plots of downstream proteins for each stress response showing various responses depending on the exposure magnitude. GFP intensity and % GFP positive cells were derived from the normalized value according to the method (Data Analysis Session). C) Time dynamic plots of Annexin V positive cells (i), PI positive cells (ii), and log fold change of cell number to the control at time point 0 (iii) as the part of the experimental readouts exhibiting the cellular outcomes. Shadow in the plots represent standard error mean (SEM). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) human signaling protein interactions (Table 2). Downstream signaling targets of investigated stress response pathways were connected to the node 'Apoptosis' in the model to reflect the degree of apoptosis ( Table 2). Additional details on the logic-ODE modeling framework can be found in the previous study [50,51].

Mathematical modelling -Optimisation
The customised version of logic modelling software suite CellNOptR was applied to convert network topology into a logic model description [44]. The R-package for modeling logic-based differential equation (logic-ODE) CNORode was applied to perform the optimisation. The optimisation was run on a clustered computer (Quad-Core Intel Xeon E5530 16 GB Memory) in 3 successive rounds for 3 h each to ensure convergence. The best results from 5 independent runs were chosen for further subsequent analyses.

Mathematical modelling -Post-optimization analysis
Optimized parameters were log10 transformed and subsequently clustered and plotted as heatmaps using the function heatmap.2 from the package ggplots with Euclidean distance in log10 scale as a distance measure. A local parameter sensitivity analysis was performed by perturbing the parameter values by 10% in both directionalities and calculating the mean of apoptotic changes across all drug concentrations being investigated.

The temporal cellular stress response activation by DILI compounds
We first mapped the long term temporal dynamics (until 60 h) of cellular stress response activation upon exposure to a set of 14 compounds known to induce DILI. In addition we included 7 compounds that do not induce DILI together with additional 8 positive reference compounds that activate specific stress response pathways. We . The different reporter cell lines were first treated with a concentration range of all the selected reference compounds to assess the feasibility to determine GFP induction over pro-longed time (Fig. 1A). Automated  image analysis pipelines allowed the quantification of total GFP intensity at the single cell level in the nucleus for the NRF2 (NFE2L2), ATF4, CHOP (DDIT3), XBP1, TP53, P21 (CDKN1A), and MDM2 reporter cell lines and in the cytoplasm for the SRXN1, HMOX1, BIP (HSPA5), BTG2, ICAM1, and A20 (TNFAIP3) reporter cell lines. TNFα was added 8 h post-exposure to the IR reporter cell lines (A20 and ICAM1) in order to induce a NF-κB-mediated inflammatory response. Each cellular stress response reporter exhibited different temporal-and concentrationdependent response patterns for the model compounds (Fig. 1B). Treatment with diethyl maleate (DEM) increased expression of SRXN1-GFP within the first 20 h, followed by a stable protein level in the remaining 40 h of imaging. Consistent with the previous study [54], the response of CHOP-GFP, a well-known effector molecule of the UPR response, reached a peak between 10 and 20 h after tunicamycin exposure followed by a decrease of the protein level depending on the concentration of the compound. The dynamics of P21-GFP induction by etoposide exhibited a rapid increase in the first 20 h and remained high during the entire exposure period. The onset and amount of expression of the different reporter proteins correlated positively with the compound concentrations. Conceptually, IR reporter cell lines were used to monitor whether xenobiotic treatment would impact on TNFα-mediated temporal dynamics of A20-GFP protein expression. Diclofenac inhibits TNFα-mediated NF-κB activation [13,37,50]. We observed that increasing diclofenac concentrations (up to approximately 400 μM) synergized with TNFα and enhanced A20-GFP levels, while very high concentrations of diclofenac reduced the expression of A20-GFP. Besides the measured average cellular GFP intensity, the percentage of cells that was committed to the cellular stress response (% GFP positive cells) was also quantified. The temporal dynamics of GFP positive cells showed a similar pattern as protein expression dynamic measured by the GFP intensity (Fig. 1B, right panel).
Our live cell imaging of the reporter cell line panels allowed a comprehensive temporal mapping of stress response activation. This mechanistic assessment was complemented with the quantification of the adverse cellular outcomes which was measured as the fraction of cell death using both propidium iodide (PI) and annexin V (AnV) as markers of necrosis and apoptosis, respectively. Most prominent cell death was observed in one of the DILI compounds, nefazodone (Fig. 1C i and ii). The onset and rate of cell death was both time and concentration dependent. At high concentrations of nefazodone (>150 µM), approximately 50% of the cells died within 24 h, while at lower concentrations (less than150 µM), less than 50% of the cells died at later time points (up to 60 h). As another indicator of cell health, based on nuclear Hoechst 33342 staining, the cell number was also quantified over time (Fig. 1C  iii). At high concentrations, the cell number decreased over time, proportionally to the increase of the cell death, consistent with the decrease in cell number caused by cell death. Overall, long term live imaging of our reporter cell lines enabled us to delineate the temporal dynamics of stress response pathways as well as the corresponding adverse cellular cytotoxicity outcomes when exposed to xenobiotic.

Comprehensive mapping of DILI compound based on temporal cellular stress response activation.
Next, we quantified the effect of all 30 compounds (DILI, non-DILI, positive references; Table 1) on the activity of all 13 GFP cellular stress reporters. To establish simplified concentration temporal response maps for comparison, the 60-hour temporal dynamics (response) based on GFP intensity, fraction of AnV and PI positive cells, and cell counts (features) for each reporter cell line were transformed into an individual map with blue to red color depicting the extend of the responses for every single compound ( Fig. 2A). The color intensity indicates the magnitude of the response whether induced or positive (red) or repressed and negative (blue). A comprehensive heatmap displaying the overview of the stress responses and adverse cellular outcomes for each compound exposure was thereafter created from each individual map (Fig. 2B). Based on all the individual compound stress response and adverse cellular outcomes maps, a hierarchical clustering analysis was performed from which we identified five sub-groups of compounds activating a set of reporters in a similar fashion. Within these groups, we also identified set of compounds that also induce cell death classified in area 6 (Fig. 2B). Fig. 2C is a simplification of the clustering where the main branches in the tree structure show how certain DILI compounds cluster together with the model compounds based on specific reporter activation (branch 1: etoposide; branch 2: tunicamycin; branch 4: DEM). Interestingly, most negative controls and solvent controls are clustered together in the third area/branch. Also cyclosporine A and acetaminophen cluster in this third branch, likely due to limited overall impact on cellular stress response activation and/or limited bioactivation to toxic reactive metabolites and cytotoxicity onset.
As expected based on the clustering with the model compounds, the different branches in the dendrogram were associated with specific stress response clusters of the heatmap (Fig. 2B): branch 1 with DDR activity; branch 2 with UPR activity, and for several compounds (nefazodone, troglitazone, nitrofurantoin, and CDDO) also OSR activity; branch 3 without substantial activity in any of the reporters; branch 4 with OSR activity alone; branch 5 with low-to-intermediate activity of multiple cellular stress response pathways and limited effect on cellular outcome. Compounds within clusters 1 and 2 clearly led to adverse cellular fates (cell death induction associated with a reduced cell count). Modulation of IR activity was not a strong determinant for the clustering of DILI compounds; indeed, diclofenac, the reference DILI compound for modulation of TNFα-induced IR clustered together with low-tointermediate stress response-inducing compounds in the fifth cluster. This was associated with activation of other stress pathways by diclofenac. We also find consistent observations from the results derived from the fraction GFP positive cells (data not shown).

UPR-and OSR-related reporter activities correlate with the onset of cell death
As a next step we aimed to systematically determine the correlation between reporter activities with cytotoxic responses. In order to reduce the dimensionality of the dataset, the dynamic response of every cell line was summarized as one unique data point: the area under the curve (AUC) value. The AUC was calculated for each compound, each concentration, each reporter cell line and cellular outcomes. Thereafter, a correlation analysis was performed to determine the correlation coefficient between all parameters (cellular responses) in the matrix (Fig. 3A). As expected, the correlation analysis showed high correlation between reporters of the same stress response pathway (e.g. BTG2-GFP and TP53-GFP (R2 = 0.94); CHOP-GFP and ATF4-GFP (R2 = 0.85)). Also, AnV and PI showed a high correlation (R2 = 0.93) indicating concordance between both readouts. Cross-pathway correlation analysis revealed interesting relationships. For example, we also observed that NRF2-GFP levels showed higher correlations with ATF4-GFP (R2 = 0.64) and CHOP-GFP (R2 = 0.55) than with SRXN1-GFP (R = 0.35), despite the fact that SRNX1-GFP expression is dependent on NRF2 activity [48]. Across all pathways, ATF4-GFP, CHOP-GFP and NRF2-GFP reporter activity had the highest correlation with cell death (AnV and PI) with the most negative correlation with the cell count. We further focused on the correlations of these three proteins with the cell death outcomes by plotting all individual AUC values (Fig. 3B). The high correlation coefficients between ATF4, CHOP, and NRF2 with cell death features were mainly caused by a selected set of compounds including nitrofurantoin, nefazodone, and troglitazone at middle-to-high concentrations. These compounds were also shown to induce UPR and OSR with high cell death responses as demonstrated in Fig. 2B. Altogether, these findings suggest that there is a high correlation between the activation of UPR and OSR stress response pathways and the cell death responses.

Logic-ODE modeling captures relationships between UPR, OSR and cell death outcomes
To investigate the signal transduction of stress response pathways which governs cell death upon toxic compound treatments, we applied a mechanistic modeling analysis using the Logic-ODE approach. First, we built a minimal topological model based on the canonical structure of the four stress response pathways with crosstalk interactions previously described in a curated human signaling network compendium 'Omnipath' (see Methods). We subsequently fitted the model to normalized experimental data from the time-resolved dose-response GFP reporters where all parameters were optimized for each individual compound (see Material and Methods). Model structure, model fit quality, and optimized parameters of the Logic-ODE models are shown in Fig. 4. The majority of the Logic-ODE models (17 out of 21 compounds - Fig. 4B) fitted reasonably well to our current model topology Fig. 4A with fitting costs represented by mean square error (MSE) of less than 0.003. The model of nefazodone for instance had the highest fitting cost (close to 0.01), but still fitted relatively well to the complete GFP dataset (data now shown).
As a next step, we proceeded with subsequent parameter analysis and in silico prediction from the fitted model. In the framework of logic-ODE, an interaction between two molecules is defined as an 'edge' where each molecule/protein in the model is referred as a 'node'. The higher the parameter values for edges and nodes, the more influential the respective connections and molecules become. According to the model optimization results, the clustering of the edge parameters 'k' (Fig. 4C), which represent the speed of interactions, and the node parameters 'tau', which represent the responsiveness of a molecule to an activation (Fig. 4D), both showed that the DILI compounds are mostly clustered together (dendrograms not shown). Ascorbic acid (ASC), buspirone (BUS) and verapamil (VPL), while labelled as negative controls for DILI, still activate stress response pathways and apoptosis, hence their parameters clustered together with the ones of DILI. Among the 'k' parameters (Fig. 4C), two main clusters are observed (see dendrograms). Parameters on the cluster to the right, which comprise the interaction linking CHOP and NRF2 towards apoptosis, are active in almost all models of DILI compounds, highlighting the potential roles of UPR and OSR pathways in modulating the onset and/or predicting the liability of apoptosis. Similarly, the clustering of 'tau' parameters ( Fig. 4D) demonstrates that the node 'Apoptosis' is clustered together with the nodes from the ATF4-CHOP axis of the UPR as well as to the OSR nodes. These findings underline the potential regulatory role of these two stress response pathways in determining cell fate upon treatment with DILI compounds.

ATF4-CHOP and NRF2 branches are predicted to be main regulators of cell death outcomes
Next, we performed a sensitivity analysis to evaluate the effect of node and edge parameters on the apoptosis outcome once their parameter values are varied (e.g. a higher edge parameter 'k' implies a reaction rate increase). A local parameter sensitivity analysis was performed by perturbing the optimized parameter values by 10%, and we monitored the mean of percent apoptosis changes in response to these in silico perturbations. Fig. 5A shows top 3 parameters which are most sensitive towards the modulation of apoptotic level in each compound and we found that they are mostly the components of the UPR pathway, especially the ATF4-CHOP-Apoptosis axis. Other parameters having a strong regulatory role in apoptosis include components of the OSR pathway, especially NRF2, and also via crosstalk interactions between stress response pathways in several compounds such as diclofenac (DIC) and omeprazole (OMZ).
Among these top sensitive parameters identified from sensitivity analyses, the parameter ATF4_k_DDIT3 was the most sensitive one for nitrofurantoin (NIT). We displayed the predicted time-course dynamic of apoptosis onset after perturbing the optimized parameter values by 10% (upwards and downwards) for each NIT concentration (Fig. 5B). The effect of the perturbation on apoptosis is strongest at higher doses and at later time points. Also, the predicted apoptotic profiles do not always exhibit proportional alteration of responses after perturbing parameter values the same range. Here, decreasing of parameter values has a more dominant effect on apoptosis at high doses, potentially due to the saturation of apoptotic signal (see lines of triangles pointing downwards in Fig. 5B). Altogether, these analyses from the Logic-ODE model suggest that the UPR was the most crucial stress response pathway that controls the apoptosis in response to the investigated DILI compounds, followed by NRF2 activation. Altogehter, the output from our unbiased modelling approach is consistent with the correlation analysis based on AUC values (Fig. 3).

Functional assessment of UPR components in the control of cytotoxicity induced by DILI compounds
Finally, given the strongest predictivity of UPR pathway activation and apoptosis, we aimed to determine which UPR components are critical determinants for the onset of cell death by DILI compounds. For this, we performed a set of experiment of RNA interference-based depletion of several key UPR components (EIF2AK3, ATF4, DDIT3, XBP1, HSPA5, and ATF6). We selected compounds that induced UPR cellular stress reporter activities clustered in three different groups based on their temporal dynamics and pattern of UPR stress reporter activation in association with various degrees of cellular outcomes. For these experiments we selected one concentration per compound based on our initial screening data and that ranged between 20 and 100 Cmax, depending on the potency for induction of cytotoxicity. The CHOP-GFP HepG2 cell line was utilized for evaluating the effect of the knock down of the above mentioned UPR regulators (Fig. 6A). We use the CHOP-GFP HepG2 cell line because this protein shows the clearest response among other UPR reporter cell lines. The knock down efficiency of every different UPR component was also determined at the mRNA level and showed to be higher than 70% (Fig. 6B). Firstly, we focused on tunicamycin and cyclosporin A, as both compounds cause a rapid activation of all the UPR reporters which is sustained for prolonged time at high concentrations (Fig. 7A). siCHOP fully inhibited the CHOP-GFP induction by both tunicamycin and cyclosporin A, as expected, which did not affect the onset of cell death ( Fig. 7B and Cdark green). On the contrary, siHSPA5 caused an enhanced CHOP-GFP induction by tunicamycin and cyclosporin A (Fig. 7B -dark blue), which was associated with increased cell death in the case of TUN (Fig. 7C -dark blue). Both siELF2AK3 and siATF6 treatment enhanced cell death by tunicamycin and cyclosporin A (Fig. 7C, light green and light blue), which was also associated with differential modulation of CHOP-GFP induction dynamics compared to siControl (Fig. 7B, light green and light blue). Minuscule differences on the cell death responses of the siRNA knocked down were noticed in the DMSO treatment (data not shown).
Next, we determined the effect of depletion of the individual UPR target genes on DILI compounds that demonstrated a sustained and strong concentration-dependent induction of ATF4-GFP and CHOP-GFP without impacting on XBP1-GFP and BIP-GFP expression, including nitrofurantoin, diclofenac, nefazodone, and troglitazone (Fig. 8A). Troglitazone caused activation of all four UPR reporters, suggesting stronger cell stress responses (Fig. 8A). Depletion of in particular ATF6 sensitized cells towards nitrofurantoin-and diclofenac-induced cell death (Fig. 8B), without particularly impacting on the induction of CHOP-GFP (Fig. 8C). Nefazodone already caused quite some induction of cell death under siControl conditions, yet a trend was observed that siATF6 also promoted cell death by nefazodone (Fig. 8B). On the contrary, troglitazone also caused substantial cytotoxicity under control conditions (Fig. 8B), which was not impacted at all by depletion of any UPR modulator, suggesting critical other cellular effects that drive cell killing by troglitazone, such as likely mitochondrial damage [36,39,40].
Finally, we focused on omeprazole and azathioprine (both 100X Fig. 7. Role of individual UPR components in progressive cellular outcomes caused by tunicamycin and cyclosporin A. A) Temporal dynamics plot of the UPR components upon the compound exposures expressed by the (min-max) normalized integrated GFP intensity. Shadow in the plots represent standard error mean (SEM). B) HepG2 CHOP-GFP cells with indicated depletion of UPR components were exposed to either 10 µM tunicamycin (left) and 20 µM cyclosporine A (right). Area under the curve of cell death fraction measured from 0 to 60 h in the most affecting knocked down condition (siATF6, siELF2AK3, and siHSPA5). The error bars show standard deviation from 3 to 4 replicates. Significance is depicted with * = p-val ≤ 0.05; ** = p-val ≤ 0.01; *** = p-val ≤ 0.001. C) Temporal dynamics of CHOP-GFP upon the exposure of compounds in the UPR perturbed cells showing increase and sustained of CHOP expression (siATF6, siEIF2AK3, and siHSPA5). Shadow in the plots represents standard error mean (SEM).
(caption on next page) L.S. Wijaya et al.
Cmax) which both induced a slightly delayed onset of the activation of ATF4-GFP and CHOP-GFP at around 10 h, with omeprazole demonstrating the strongest response. Also some activation of XBP1-GFP (omeprazole only) and BIP-GFP (azathioprine only) was observed (Fig. 9A). Only for omeprazole depletion of ATF4 and ATF6 sensitized cells to enhanced cell death in association with stronger induction of CHOP-GFP; for azathioprine depletion of UPR modulators did not stimulate onset of cell death ( Fig. 9B and C).
Overall these data indicate differential temporal dynamics and selective activation of different UPR reporters by compounds with DILI liabilities. Moreover, in particular ATF6 stands out as an important adaptive modulator to repress onset of cell death for various DILI compounds that strongly activate the ATF4-CHOP axis of the UPR.

Discussion
In this study, we have demonstrated the feasibility to integrate the temporal dynamics of stress response pathway activation and the correlation with the progressive cellular outcomes utilizing combined experimental and computational approaches. Here, we have established comprehensive maps of the temporal dynamic activation of cellular stress response pathways including OSR, UPR, DDR, and IR as well as cell death information caused by various DILI compounds. These datasets were used for builindg a Logic-ODE modelling to uncover the most critical cellular stress pathways that are predictive for DILI outcome. The major finding was that the UPR on ATF4-CHOP axis is most strongly associated with DILI compound induced cytotoxicity. The role of UPR involvement in DILI-induced cytotoxicity was further established by RNA interference experiments, and identified ATF6 as the most critical determinant of DILI compound induced cytotoxicity.
Our experimental data demonstrate the different dynamics of the proteins involved in the specific stress response pathways by various DILI compounds. The different responses were found to be affected by the degree and duration of exposure. Mostly, the responses were clearly concentration dependent: both activation rate and extend of the response increased with the compound concentration. However, for some compounds such as cisplatin, nefazodone, and nitrofurantoin, the high concentrations paradoxically showed limited activation of particular stress response pathways. This is most likely due to rapid progressive adverse cellular responses that may have impacted on general transcription and translation machineries and, thereby, affected the expression of newly synthesized GFP reporter proteins. Despite this general principle of stress responses activation, each compound elicited distinct temporal dynamics. Clustering based on those response dynamics and cellular outcomes showed that the DILI compounds activating UPR induced strong progressive cell death responses. Meanwhile, the compounds inducing OSR without the co-activation of other stress responses did not lead to the adverse cellular outcome such as the non-DILI control compounds. Interestingly, some compounds showing progressive outcomes induced both UPR and OSR, particularly NRF2 nuclear translocation. The activation of OSR by these compounds might be involved in the partial protection of the cells. Of relevance, we found moderate correlation between NRF2 and SRXN1 as its downstream target. This low correlation is linked to the different temporal dynamics of these 2 proteins: the nuclear translocation of NRF2 typically occurs early within the first period of compound exposure followed by a decrease of this protein while SRXN1 exhibits a rather sustained induction over time with a plateau level after twenty hours depending on the compound. Furthermore, we found that there is a positive correlation between the activation of UPR components ATF4 and CHOP with cell death. This is in line with previous studies which showed that UPRactivating compounds are mostly classified as the severe DILI compounds [7,13].
In order to get a deeper mechanistic insight of our experimental findings, we applied a Logic-ODE approach to describe the mode of stress response pathways regulation upon DILI compound perturbation. Unlike traditional mechanistic ODE-based models, this modeling approach requires only prior knowledge on signaling network architecture, without information from biochemical reactions or kinetic parameters which are often unavailable. In this work, we demonstrate that the Logic-ODE approach allowed us to capture changes in molecular responses in both concentration-and time-response fashion from our dynamic experimental data which is similar to what traditional ODEbased models can deliver [37,54]. With the simplicity of the Logic-ODE framework, it was also possible to build customized models for all investigated compounds.
Apart from linking the parameters from Logic-ODE models as features to predict pathological observations [12], we are the first to demonstrate that the model states in the Logic-ODE models could also be directly linked to cell phenotype such as apoptosis. This novel implementation allows to mechanistically connect the modulation of stress response signals towards cytotoxicity measures. In addition, standard modeling analytical techniques such as parameter sensitivity analyses could be performed to predict the changes of molecular responses due to the modification of model parameters, allowing the sensitive parts of the network that control progressive cellular outcomes to be identified. Furthermore, as the expected profiles of signaling molecules and cellular outcomes can be predicted over the range of concentrations and timecourse, the Logic-ODE models could further guide the selection of compound doses and time-points for a validation experiment which demonstrates the highest effects from the respective perturbation that can still be observed experimentally.
Our unbiased Logic-ODE modelling identified the parameters from the ATF4-CHOP axis of the UPR pathway with the highest correlation with cell death/apoptosis outcome and highest sensitivity across multiple DILI compounds (Fig. 5A). These modeling results were in agreement with the correlation plots of individual reporters with adverse outcomes (Fig. 3), providing confidence in the Logic-ODE modelling approach. A role for the ATF4-CHOP axis in the control of apoptosis is in agreement with the literature [11,15,30,33]. Analysis of a large DILI compound transcriptomics dataset in primary human hepatocytes has confirmed that various DILI compounds induce the activation of the ATF4-CHOP axis [13]. However, a systematic assessment on the involvement of UPR components in the control of DILI compoundinduced cytotoxicity has to our knowledge never been performed. Neither CHOP nor ATF4 were critical determinants of DILI compoundinduced cytotoxic outcome, possibly due to other compensating cytoprotective pathways arising upon the depletion of these proteins such as increase of receptor-mediated signaling pathways (Niemeijer et al. unpublished observations). However, we identified ATF6 as a critical determinant of the onset of apoptosis by several DILI compounds that also demonstrated a significant activation of the UPR response, Fig. 8. Role of individual UPR components in progressive cellular outcomes caused by nitrofurantoin, diclofenac, nefazodone and troglitazone. A) Temporal dynamics plot of the UPR components upon the compound exposures expressed by the (min-max) normalized integrated GFP intensity. Shadow in the plots represent standard error mean (SEM). B) HepG2 CHOP-GFP cells with depletion of indicated UPR components were exposed to 234 µM nitrofurantoin, 808 µM diclofenac, 79 µM nefazodone and 381 µM troglitazone. Area under the curve of cell death fraction measured from 0 to 60 h in the most affecting knock down condition (siATF4 and siATF6). The error bars show standard deviation from 3 to 4 replicates. Significance is depicted with * = p-val ≤ 0.05; ** = p-val ≤ 0.01; *** = p-val ≤ 0.001. C) Temporal dynamics of CHOP-GFP (expressed by the (min-max) normalized integrated GFP intensity) upon the exposure of compounds in the UPR perturbed cells increase and sustained of CHOP expression (siATF4 and siATF6). Shadow in the plots represents standard error mean (SEM).
(caption on next page) L.S. Wijaya et al. including cyclosporin A, nitrofurantoin, diclofenac and omeprazole. Yet, cytotoxicity caused by other DILI compounds, troglitazone, nefazodone and azathioprine, that also demonstrated activation of the UPR pathway reporters, was not significantly modulated by depletion of any UPR component; although nefazodone showed some trend that mimicked the responses from nitrofurantoin and diclofenac. This suggests that despite onset of UPR by these latter compounds, that other cellular injuries, e.g. mitochondrial toxicity, may dominate the adverse cell fate. Recently, activation of the ATF4 pathway by potent mitochondrial respiratory chain inhibitors was demonstrated [28], suggesting that a contribution of mitotoxicity by some of the DILI compounds in the activation of the ATF4-CHOP pathway cannot be excluded. Overall, these findings highlight the fact that there are still missing regulatory mechanisms within the UPR branches which could not be captured in this particular study.
We observed that in particular, ATF6 depletion resulted in a significant increase of cell death in the majority of the tested DILI compounds. This might imply an adaptive protective role of ATF6 in DILI compoundinduced UPR activation. Moreover, activation of ATF6 is known to increased BIP expression (where silencing BIP expression in this study also increases cell sensitivity) thereby providing renewed cellular proteostasis and cytoprotection [3]. This is in line with previous studies reporting that blocking the ATF6 pathway leads to the increase of cell sensitivity towards apoptosis during stress [43]. Previously, we established a key role for ATF6 in modulating the CHOP response of tunicamycin [54]. In this study we showed that depletion of ATF6 impacts considerably on the overall CHOP-GFP expression upon DILI compound exposures by increasing and/or stabilizing its expression. This might indicate that loss of ATF6 increase the sensitivity of the cells by prolonging the ATF4-CHOP axis activation. This finding is supported by a previous study showing that chronic UPR activation indicated by CHOP expression induces cytotoxicity [17]. Despite the concordance to the previous studies, the usage HepG2 cell lines as the most high throughput liver model in this study involves some limitation. The difference biological properties from primary cells such as proliferation rates, central metabolism and biotransformation capacity, could lead to different response compared to PHH or in vivo situation. This might be the reason of low-mild activation of stress responses in negative control compounds such as verapamil. The low metabolic capacity of HepG2 cells lines might also arise the over/underestimation of toxic responses resulted by the toxic metabolites as observed for acetaminophen. These issues could be overcome by culturing the HepG2 in spheroids as it was previously reported that in the 3D spheroid culture with higher CYP450 [22] or cultured in medium for improved differentiation [4].
In conclusion, we presented a new and comprehensive analytical approach by combining highly granular dynamic experimental data of cellular stress response pathway activation with Logic-ODE modeling. Our integrated approach indicates that activation of UPR on ATF4-CHOP axis is an important predictor for adverse cellular outcomes for DILI compounds. We propose that our panel of UPR cellular stress response reporters and the temporal dynamic analysis of the UPR on ATF4-CHOP axis activation can be a valuable component of a DILI derisking strategy for novel drug candidates. Fig. 9. Role of individual UPR components in progressive cellular outcomes caused by omeprazole and azathioprine. A) Temporal dynamics plot of the UPR components upon the compound exposures expressed by the (min-max) normalized integrated GFP intensity. Shadow in the plots represent standard error mean (SEM). B) HepG2 CHOP-GFP cells with depletion of different indicated UPR regulators were exposed to 470 µM omeprazole (left) or 34 µM azathioprine (left). Area under the curve of cell death fraction measured from 0 to 60 h in the most affecting knocked down condition (siATF4 and siATF6). The error bars show standard deviation from 3 to 4 replicates. Significance is depicted with * = p-val ≤ 0.05; ** = p-val ≤ 0.01; *** = p-val ≤ 0.001. C) Temporal dynamics of CHOP-GFP (expressed by the (min-max) normalized integrated GFP intensity) upon the exposure of compounds in the UPR perturbed cells showing increase and sustained of CHOP expression (siATF4 and siATF6). Shadow in the plots represents standard error mean (SEM).