Unraveling cellular pathways contributing to drug-induced liver injury by dynamical modeling

ABSTRACT

Introduction: Drug-induced liver injury (DILI) is a significant threat to human health and a major problem in drug development. It is hard to predict due to its idiosyncratic nature and which does not show up in animal trials. Hepatic adaptive stress response pathway activation is generally observed in drug-induced liver injury. Dynamical pathway modeling has the potential to foresee adverse effects of drugs before they go in trial. Ordinary differential equation modeling can offer mechanistic insight, and allows us to study the dynamical behavior of stress pathways involved in DILI.

Areas covered: This review provides an overview on the progress of the dynamical modeling of stress and death pathways pertinent to DILI, i.e. pathways relevant for oxidative stress, inflammatory stress, DNA damage, unfolded proteins, heat shock and apoptosis. We also discuss the required steps for applying such modeling to the liver.

Expert opinion: Despite the strong progress made since the turn of the century, models of stress pathways have only rarely been specifically applied to describe pathway dynamics for DILI. We argue that with minor changes, in some cases only to parameter values, many of these models can be repurposed for application in DILI research. Combining both dynamical models with in vitro testing might offer novel screening methods for the harmful side-effects of drugs.

KEYWORDS:

• Drug induced liver injury (DILI)
• dynamical modeling
• ordinary differential equations
• stress pathway
• death pathway
• predictive toxicology
• cell fate
• adverse outcome pathway

1. Introduction

Drug-induced liver injury (DILI) is the most frequent cause of acute liver failure [1]. Almost half of the cases of such liver failure in the United Kingdom and United States are attributed to DILI [1,2]. This problem represents a huge financial burden to the pharmaceutical industry, because new medication is often found to be causing DILI at a very late drug development stage, such as when a drug is already in its clinical test phase or, even worse, when the drug is already on the market. DILI is the most common reason for drug withdrawal [3] and can be separated into intrinsic or idiosyncratic liver injury. Intrinsic liver injury is dose dependent and predictable when given in sufficiently high dose. The mechanism behind intrinsic DILI frequently involves exposure of cells to an excess of stress so that a cell either dies after activating highly regulated cell death pathways (i.e. apoptosis) or directly dies because it cannot deal with the high exposure (i.e. necrosis). In contrast, idiosyncratic liver injury is unpredictable and has a low incidence rate (1:100,000 users) [4]. The mechanisms causing idiosyncratic DILI are highly complex and poorly understood. Current techniques that include animal testing or simple in vitro tests are insufficient to predict the idiosyncratic DILI properties of drugs, although many efforts to enhance in vitro testing have been started [5] (e.g. the IMI MIP-DILI project). In 2007, the National Research Council Committee of Toxicity and Assessment of Environmental Agents presented a new vision and strategic plan for toxicity testing [6]. In short, the suggested approach involves a combination of dynamic in vitro data and computational modeling.

Many types of models are of interest for the understanding and prediction of DILI, including quantitative adverse outcome pathways (qAOPs), cheminformatic models, pharmacokinetic–pharmacodynamic (PK–PD) modeling, and dynamical pathway modeling with ordinary differential equation (ODE) models. Excellent reviews have recently been published about adverse outcome pathways [7], multi-scale computational models [8], and cheminformatic models [9,10].

AOPs describe the path from a molecular initiation event to an adverse outcome, through a number of key events. These key events can occur on different biological levels (subcellular, cellular, tissue, organ, individual, and population level) (Figure 1, bottom with a microscope image reproduced from [11]). Adverse outcomes in DILI include clinical observations such as liver fibrosis, liver steatosis, and cholestatis. The main goal of constructing an AOP is to apply it in risk assessment based on mechanistic reasoning [12]. qAOPs are closely related to multi-scale computational models, which also try to bridge different scales yet often explicitly include spatial effects as well [8].

Figure 1. Overview of different scales relevant for DILI and the contribution of dynamical pathway modeling. Upon exposure to drugs, cellular stress pathways are activated to maintain homeostasis. For illustration, we show two pathways for the unfolded protein response and oxidative stress response. Malfunctioning of the pathways could result in an imbalance of cell fates. Such malfunctioning can progressively change states of the liver (or other relevant organs), eventually leading to various adverse outcomes in a clinical setting, like cholestasis. Here, dynamical modeling aims to quantitatively describe the underlying intra-cellular signaling networks. This computational approach potentially offers a mechanism-based tool towards predictive toxicology to understand and eventually to reduce the occurrence of DILI. For illustration purposes, we include a 2D representation of cells in a single microscopy image reproduced from [11] under the terms of the Creative Commons Attribution License.

Cheminformatic models, i.e. in silico models based on chemical structure, employ existing data on chemical properties, their key chemical substructures, effect, and interaction capabilities in order to make predictions for untested chemicals. They are helpful in the classification of potentially useful chemicals, or give insight in potentially harmful side effects like DILI [13].

PK–PD modeling aims to predict exposure–response relationships. Exposure is measured as the concentration of a drug in the plasma, and response is the effect the drug has on the patient. Thus, the efficiency with which the drug reaches the target location and the effect on that location are included in these models [14].

Here, we will review the recent progress in dynamical modeling of toxicity pathways, i.e. cellular response pathways that, when sufficiently perturbed, are expected to contribute to adverse health effects such as liver injury [6]. These pathways can be divided into two groups, i.e. cellular stress and death pathways. Death pathways are networks involved in the process of highly regulated cell death, whereas stress pathways are networks capable of responding quickly to a threat, resulting in the production of proteins that help to remove the threat or minimize damage. Dynamical modeling aims to describe these cellular processes as a network of biochemical reactions, thus representing the biological processes schematically and quantitatively. The biochemical reactions are formulated as ODEs such that they obey kinetic principles like the law of mass action, Michaelis–Menten kinetics, or cooperativity modeled by a Hill equation [15,16]. Such models can offer mechanistic insights about the molecular processes within a pathway. For example, the ODE models can be used to investigate the potential stable states of the system and to study what happens when a pathway is perturbed. If the system in that case evolves toward a stable malfunctioning state, this may be relevant for DILI-related disease progression. Using ODE models, one can furthermore test hypotheses about the network structure and how that structure affects the resulting dynamics of a cellular system [17]. The expected impact of specific biochemical interactions on the dynamics can be evaluated within models, which would be much harder or even impossible by experimental work. This is because the mechanisms involved are often highly nonlinear and include multiple feedback loops, which makes interpretation of experimental data complicated and sometimes counterintuitive. The aim of this review is to assess recent progress in mechanistic dynamical modeling of pathways contributing to DILI and highlight its biological relevance.

2. Stress and death pathways

Technological developments in molecular biology, e.g. next-generation sequencing and high-throughput measurements, have enabled significant advances in the field of systems biology. For example, it has resulted in new insights into formerly unknown key network elements. Indeed, the application of cluster analysis on ‘big data’ has given a better understanding of the complex network of the cell. However, this large amount of information can be daunting, and deciding which pathways and players one should focus on becomes an increasingly difficult question. This is because cluster analysis provides an overview of which players might be involved, but it does not directly provide information about their causal relations and whether direct or indirect interactions occur. Combined experimental and computational approaches are required to solve this problem [18,19].

In this review, we focus on the recent modeling progress with respect to cellular stress and death pathways that are suspected to be related to DILI, aiming to obtain a mechanistic understanding of the functioning of the individual pathways. Stress pathways try to repair a ‘controlled variable’ (e.g. temperature, levels of reactive oxygen species [ROS], or a stressor itself) that is out of control [20]. To this purpose, a stress pathway typically consists of two major components (Figure 1, top right). First, a ‘sensor’ that picks up the stress-induced deviation of the controlled variable, leading to activation of the transcription factor. Second, a transcription factor that can rapidly respond, usually because under homeostatic conditions, it is already abundantly present, albeit in an inactive state. Death pathways are often instigated by a stress pathway, so if the cell is unable to cope with the threat, it will activate a death pathway to minimize damage to surrounding cells. This avoids a spillage of the content of dead cells over neighboring cells, which could otherwise cause a harmful inflammatory response [21].

ODE-based models of cellular stress pathways contain system variables and system parameters: System variables represent concentrations of sensors, transcription factors, ‘controlled variables,’ and possibly intermediate molecules. The system parameters characterize molecular interactions and processes like degradation quantitatively. Although parameters are constant in a particular set of ODEs, their values could differ amongst cell types or due to treatment. A general mathematical formulation of ODEs is shown below the simplified control diagram of a stress pathway in Figure 1. System variables from a stress pathway can act as inputs to determine the fate of a cell, leading to an integrated model that predicts the effect of various drug treatments on cell fate.

Determining which stress and death pathways are involved in DILI is an important part of predicting potentially harmful compounds. The pathways discussed below, i.e. the apoptosis pathway and the pathways dealing with oxidative stress, inflammatory stress, DNA damage, unfolded proteins, heat shock, are frequently implicated in DILI (for a schematic overview of the most important components of these pathways, see Figure 2). However, this list is unlikely to be complete, and the search for important indicator pathways continues.

Figure 2. Schematic overview of the regulatory network underlying various stress pathways. To illustrate their adaptive anti-stress roles, only key players in the anti-stress feedback of each pathway are presented in each panel. A: Oxidative stress, B: Inflammatory stress, C: Unfolded-protein stress, D: Heat-shock stress, and E: DNA-damage stress. See the main text for details and other relevant biochemical interactions.

Apart from the search for DILI-related pathways, we advocate applying stress and death pathway modeling to liver cells in order to help understand and predict DILI. The oxidative stress and the heat-shock response (HSR) pathways have already been applied to the liver [22,23], but to our knowledge, this is not yet the case for other stress pathways. As we will discuss below in more detail, general models or models for different cell types or organs can be repurposed for liver cells. In case of the inflammatory stress pathway, this is most complicated, because collaboration of multiple cell types takes place during an inflammatory response in the liver. Still, dynamical models can give valuable mechanistic insight even if based on only one cell type. Linking stress pathways to death pathways can be an important step toward predicting the severity of DILI. Such linking of pathways has been performed in several modeling studies (discussed in more detail below) but has not yet been applied to liver cells [24,25].

2.1. Oxidative stress pathway

The oxidative stress pathway is considered a highly important pathway for protection against drugs that can lead to liver injury [26]. This pathway responds to ROS and/or electrophiles. ROS are reactive molecules containing oxygen, e.g. reactive drug metabolites or hydrogen peroxide produced by mitochondria. Under homeostatic conditions, mitochondria produce a limited amount of ROS, but under stressful conditions, the ROS production rate increases. In response, the oxidative stress response pathway promotes the expression of proteins that deactivate ROS, thus preventing further cell damage caused by ROS. When ROS are not adequately neutralized, other stress pathways likely also get activated because the sustained high levels of ROS leads to DNA damage as well as build-up of unfolded proteins due to oxidation.

2.1.1. Mechanism

Nuclear factor erythroid 2-related factor 2 (Nrf2) and Kelch-like erythroid cell-derived protein with CNC homology-associated protein 1 (Keap1) are the two major players in this stress pathway (Figure 2(a)). Nrf2 is a transcription factor that regulates genes involved in the neutralization of oxidative stress. Keap1, an inhibitor of Nrf2, is very sensitive to ROS and represents the ‘sensor’ of the pathway. ROS and electrophiles are the ‘controlled variables’ of the system, they trigger the pathway when their concentration exceeds a threshold, and are managed by Nrf2-regulated genes. Recently, a theory has been proposed regarding the mechanism explaining the functioning of the Nrf2–Keap1 complex called the Hinge and Latch theory [27]. This theory proposes that in the absence of oxidative stress, Keap1 induces Nrf2 degradation, which is compensated for by a high Nrf2 production rate. When Keap1 interacts with ROS, it changes its binding to Nrf2 and as a consequence no longer promotes degradation. Therefore, Nrf2 proteins increase in numbers in the cytoplasm, move freely into the nucleus, and activate genes that help to resolve oxidative stress. Several important targets of Nrf2 are involved in glutathione (GSH) synthesis, thus ensuring increased levels of cellular GSH and hence increased neutralizing capacity for electrophiles and ROS.

2.1.2. Modeling

Although the Nrf2 pathway is widely recognized as an important stress pathway, only few dynamical models have been developed for it. The first model describing the oxidative stress pathway was constructed by Zhang et al. [20], using a homeostatic control system as a basis. This model describes in detail how the Nrf2 pathway responds to stress and how it alters gene expression levels, and ultimately deals with ROS via GSH activity. The general aim was to describe antistress gene regulatory networks, to explore the ‘design principles’ that have evolved, and use the oxidative stress pathway as an example. One of the design principles that was investigated involves the response coefficient [20,28,29].

A response coefficient (also called gain) defines the sensitivity of a system variable of interest to a perturbation when the system is at equilibrium. In the case of stress pathways, the response coefficient indicates how much a stressor influences a stress-relieving molecule. A high response coefficient indicates that the stress pathway can respond fiercely to only a small amount of a stressor. There are multiple ways to enhance the response coefficient within a stress pathway. First, the homo-dimerization or homo-trimerization can increase the response coefficient, which frequently occurs for transcription factors before they become transcriptionally active or for reactions involving enzymatic complexes. For example, the oxidative stress response pathway regulates two homodimer antioxidant enzymes (GSH reductase and superoxide dismutase). Second, the response coefficient can be enhanced by autoregulation, for instance when a transcription factor induces its own production. For example, Nrf2 enhances its own production and therefore regulates itself. Stress pathways often combine multiple ways to enhance the response coefficient, which strongly affects the shape of the dose response curve [20].

Zhang et al.’s [20,28,29] findings highlight why it is in general important to apply computational modeling to stress pathways on top of experimental work. With respect to the application to oxidative stress in particular, their models were subsequently adjusted to be used for renal cells, in combination with transcriptomic, proteomic, and metabolomic data [30]. The study focused specifically on the effect that the toxicant cyclosporine A (CsA) has on renal cells. Moreover, Hamon et al. [30] extended the model with an in vitro PK model creating a pharmacokinetic systems biology (PK-SB) model [31] to take into account how much CsA actually enters the cell.

Later, the same group [22] adapted this PK-SB model for liver cells, by recalibrating the model using transcriptomic data, including ROS data at different time points and literature-based GSH data. They used the compound Acetyl-para-aminophenol (APAP) to induce cytotoxicity in a liver-on-a-chip and compared the results with the PK-SB model. The saturation of GSH metabolism led to time- and dose-dependent APAP toxicity that was in agreement with their in vitro experiments. This work demonstrates that dynamical modeling incorporating stress pathway can be used for DILI screening.

2.2. Inflammatory stress pathway

In early DILI research, drugs were classified depending on whether the immune system was involved in the response or not, dividing DILI into allergic and nonallergic responses. However, this classification turned out not to be fully accurate, because the immune response was still involved in cases that were previously classified as nonallergic. When a cell dies, especially if this is not through apoptosis, it triggers the adaptive immune response, a process called ‘sterile inflammation.’ The inflammatory stress response has been suggested to play a role in the unpredictability of DILI, and the nonlinearity of drug response curves [1,32,33].

The inflammatory response gets activated when innate immune cells detect an infection or tissue injury. To communicate to the surrounding cells that there is an infection or damaged tissue, cytokines such as tumor necrosis factor α (TNFα) is released. These pro-inflammatory cytokines are detected by receptors on both immune cells and hepatocytes. In both cell types, this causes cytokine-induced signaling. Within immune cells, this represents one of the required signals for activation of the adaptive immune response, which induces a classical allergic DILI response [34].

2.2.1. Mechanism

The main transcription factor of the pro-inflammatory response pathway is nuclear factor kappa-light-chain-enhancer of activated B cells (NF-κB) (Figure 2(b)). To be precise, NF-κB is a family of dimeric transcription factors, of which the RelA(p65):p50 complex is the most common. They regulate genes that in turn regulate a wide variety of processes such as cell division, apoptosis, and inflammation [35]. Under homeostatic conditions, the NF-κB dimers are bound to one of multiple possible NF-κB inhibitor (IκB) isoforms (the sensor), which keeps them in the cytoplasm. When the pathway becomes activated, the enzyme IκB kinase (IKK) phosphorylates IκB, which as a result becomes targeted for degradation by the proteasome. The freed NF-κB can then migrate into the nucleus and activate many genes involved in a wide variety of functions, ranging from inter- and intracellular communication, changes in cell behavior and cell fate. Some of these genes are responsible for positive (TNFα) or negative feedback loops (IκB isoforms and the A20 inhibitor of IKK).

2.2.2. Modeling

Compared to oxidative stress, the inflammatory stress pathway is much more extensively modeled and investigated. Williams et al. have already written excellent reviews about the progress on this subject [36,37]. Therefore, we will here discuss only some highlights.

NF-κB pathway modeling started at the beginning of this century by Carlotti et al. [38] who developed a simple model to quantify the shuttling speed of NF-κB and IκB in and out of the nucleus. A more complex pathway model was developed by Hoffman et al. [39] that consisted of a set of 24 ODEs. Their study focused on the regulation of three IκB isoforms (IκBα, IκBβ, and IκBε), exhibiting distinct functions. Employing a combination of mouse experiments and modeling, they showed that IκBα, through its strong negative feedback, is responsible for a rapid NF-κB response and a highly oscillatory profile, while IκBβ and ε respond more slowly and are responsible for the dampening of the oscillations. In follow-up work, Hoffman and colleagues have focused on different aspects of the pathway. For example, some of their latest work concerned NF-κB dimer formation and the role of IκBβ as a stabilizer [40], and the role of Toll/interleukin-1 receptor-domain-containing adapter-inducing interferon-β (TRIF) as a posttranscriptional regulator of TNFα [41].

The NF-κB signaling response is not homogenous throughout the population, and distinct subpopulations have been detected within the same experiment [42]. Both Tay et al. and Turner et al. have shown that not all cells respond to low concentrations of TNFα, and that this is probably due to a heterogonous digital process. IKK (the inhibitor of IκB) needs multiple phosphorylation steps to become activated, which could lead to a digital cellular response profile [43,44].

In the modeling approaches, some groups have included only parts of the NF-κB network that are sufficient to explain the dynamical behavior measured in in vitro experiments [45–48] or have developed methods to reduce complex models [49]. Krishna et al. [45] showed that three ODEs with a negative feedback loop are sufficient to simulate the observed oscillatory behavior. Zambrano et al. [48] also developed a relatively simple model with 11 ODEs and 14 parameters. They demonstrated that the observed shape of oscillations amongst different experimental conditions and cell types (e.g. spiky or smooth) could all be achieved in their model by varying multiple internal parameters of their model by maximally a factor two [48].

The NF-κB pathway is difficult to fully grasp. Many modeling studies have been performed and over and over again, the system turns out to be more complicated than expected. This is because the pathway is highly nonlinear, has multiple positive and negative feedback loops, and acts as an information hub, receiving and passing on information to many other pathways. NF-κB pathway models have to our knowledge not been applied to liver cell data. Instead, most models have been applied to immune cells, because the cytokine TNFα is important for regulation of the immune response and is itself regulated by NF-κB. For the future application of NF-κB models in DILI research, it should be kept in mind that cell types differ in the cytokines they produce and respond to. Indeed, the liver is a complex organ consisting of multiple cell types, all responding differently to cytokines such as TNFα [50].

2.3. Unfolded protein response pathway

The endoplasmic reticulum (ER) is the main organelle to fold polypeptides or misfolded proteins toward structures with their particular functions [51]. The ER is equipped with three main signaling branches to cope with stress when un- or misfolded proteins (i.e. the ‘controlled variable’) accumulate due to internal or external perturbations like aging [52] and drug treatments. Collectively, these branches are known as the unfolded protein response (UPR). In normal situations, when an overload of ER stress becomes intolerable for a cell, its intracellular network invokes the pathways for apoptotic death to prevent adverse events at tissue level [51]. Malfunctioning of the ER stress pathway can lead to prevailing apoptosis-promoting actions instead of survival, i.e. an imbalance that can also result in liver injuries [53].

2.3.1. Mechanism

Three transmembrane proteins have the ability to sense the presence of unfolded proteins: inositol-requiring enzyme 1 (IRE1), Protein kinase RNA-activated-like ER kinase (PERK), and activating transcription factor 6 (ATF6) (Figure 2C). In the absence of ER stress, each of these three proteins forms a complex with the binding immunoglobulin protein (BiP). Upon ER stress, BiP molecules dissociate from the three sensors and act as molecular chaperones of unfolded proteins, promoting their correct folding. The dissociation will also trigger regulatory mechanisms in each of the three branches: In the first UPR branch, free IRE1 molecules dimerize and auto-transphosphorylate [54]. This activation leads to alternative splicing of X-box binding protein 1 (XBP1) mRNA [55]. This mRNA activates the XBP1’s transcription factor, which is responsible for the regulation of genes involved in quality control of proteins, in the degradation of misfolded proteins, and in ER expansion. In the second branch, free PERKs first form dimers [56]. Then, PERK dimers phosphorylate the α subunit of eukaryotic initiation factor 2 (eIF2α). Phosphorylated eIF2α on the one hand attenuates general translation [57] and on the other hand activates the transcription factor ATF4. In the third branch, ATF6 first translocates to the Golgi apparatus membrane and gets cleaved [58]. This releases a cytosolic 50-kDa domain that operates as a transcription factor (i.e. cleaved ATF6 [ATF6f]) regulating quality control of proteins. The three factors, i.e. XBP1, ATF4, and ATF6f, enhance the production of BiP, which binds to newly synthesized proteins to recover the cell from ER stress [51,59]. On the other hand, the three factors also increase the translation of CCAAT/enhancer-binding protein homologous protein (CHOP). CHOP has a proapoptotic role yet also upregulates the growth arrest and DNA damage-inducible protein (GADD34). GADD34 provides negative feedback by promoting the dephosphorylation of eIF2α, which restores the protein translation capacity as soon as stress is alleviated [60]. Amongst these proapoptotic and pro-adaptive pathway branches, a quantitative understanding of the cooperation between the branches is still lacking.

2.3.2. Modeling

A previous review [61] focused on both static and dynamical models for the optimization of industrial protein production, including UPR models. Here, we review dynamical features of the UPR and existing models including reaction kinetics, cross talk between the three branches, and potential links to death pathways. A simple dynamical model of the UPR involved a linear system with separate proapoptotic (CHOP/GADD34) and pro-adaptive (BiP) variables [62]. Simulation results showed a faster transient response of the proapoptotic variable than the pro-adaptive one, which is in agreement with experimental observations. Other kinetic models of the ER stress pathway have employed nonlinear chemical reactions. First, to improve our understanding of the role of BiP binding to IRE1, the mechanism to sense ER stress by IRE1 was investigated in a detailed dynamical model [63]. Together with experimental validation, it was demonstrated that BiP provides a buffer for inactive IRE1 molecules, thus regulating IRE1 activity and protein folding homeostasis. Second, several studies investigated the IRE1 and PERK branches of the UPR and the role of attenuation of translation when there is repeated exposure to ER stress [64,65]. An example of such repeated exposures naturally occurs in secretory cells, e.g. pancreatic β-cells, in which highly variable patterns of blood glucose level cause high-frequency pulses of insulin production that are processed by the ER. Type 2 diabetic patients often exhibit irregularity of insulin pulsatility [66]. The modeling studies showed that attenuation of translation is beneficial for cells in case of repeated exposure [64], yet becomes dispensable in case of a high-frequency stress-induced stimulus [65]. With respect to the observed irregularity of insulin pulsatility, this model offered a tool to further investigate the potential role of UPR in diabetes.

More recently, computational models incorporating all three UPR branches have been constructed [24,67]. Mathematical analysis on one of these models has revealed three distinct types of dynamical behaviors of CHOP [24], i.e. transient activation, sustained oscillations, and elevated activity. Interestingly, a parameter regime with oscillations in CHOP and also in the eIF2α-mediated attenuation of translation was identified in these simulations. To date, such oscillations have not been observed in experiments, which could mean that the parameters for the oscillatory regime are not biologically realistic. These dynamical models will be useful in future exploratory studies looking into cross talk amongst the three branches.

Apart from investigating the intrinsic complexity of the biochemical network, ER stress models have been exploited to understand how different types of pharmaceutical agents induce ER stress, which can rely on different mechanisms. For example, dithiothreitol (DTT) breaks the disulfide bonds of proteins, while thapsigargin (TG) depletes Ca²⁺ from the ER [68], thus reducing the ER’s capacity to fold newly translated polypeptides. In one modeling study [65], the bond-breaking stress, e.g. induced by DTT, was considered to function in an additive way with the basal influx rate of new unfolded proteins, while the translation stress as e.g. induced by TG was considered to function in a multiplicative way with that rate. In order to study the role of translation attenuation (by phosphorylated eIF2α), two model variants were studied: one with translation attenuation upon exposure to stress and one without such attenuation. Both for translation stress (like TG) and for chemical stress (like DTT), the simulations revealed a stronger reduction of the amount of unfolded proteins in the presence than in the absence of translation attenuation. Moreover, this beneficial role of translation attenuation was most prominent for translation stress. As a final contribution of modeling studies to the understanding of the UPR, the effects of ER stressors on cell death have been investigated. Erguler et al. [24] integrated their ER stress pathway model with a previously constructed apoptosis module [69] via bcl-2-like protein 4 (also known as BAX)/Bcl-2 homologous antagonist/killer (also known as BAK1) in mitochondrial membranes. In the model, the CHOP activity represented the stress variable. This study demonstrated that a preconditioning stress at a mild level can suppress the apoptotic signal (i.e. BAX) for subsequent ER stress, since BiP molecules already accumulate sufficiently. Further modeling studies combined with quantitative data on UPR activation and cytotoxicity measurements are likely to advance our understanding of the role of ER stress in inducing or preventing cell death. Moreover, integrating such (intra)cellular models into computational models at a larger biological scale will help to better predict DILI in the future.

2.4. HSR pathway

A common environmental perturbation, heat shock, i.e. an elevated temperature, can denature proteins, potentially leading to proteotoxicity. Upon such heat shock, heat-shock proteins (HSPs) act as molecular chaperones that bind to denatured proteins (i.e. the ‘controlled variable’) and recover those into folded structures [70]. The resulting transcriptional and (post)translational processes are collectively termed the HSR (Figure 2D) and malfunctioning of the HSR can be involved in DILI [71,72].

2.4.1. Mechanism

HSPs occur in both the cytosol and in various organelles to assist protein folding (e.g. BiP in the above-mentioned ER stress pathway is also a HSP). In unstressed situations, HSPs bind to heat-shock transcription factors (HSFs), amongst which HSF1 is the major one [73]. When large numbers of denatured proteins compete with HSF1 for binding to HSP, HSF1 monomers form dimers and trimers. Subsequently, HSF1 trimers bind to the heat-shock sequence element and activate the HSP encoding gene. This leads to increased transcription of HSPs which bind to remaining denatured proteins to resolve the stress [70]. As discussed above, the multimerization of HSF1 monomers increases the gain of the response, thus strengthening the protective capability to small thermal fluctuations [20].

2.4.2. Modeling

Dynamical models of the HSR have been built and studied extensively with data sets in hyperthermia conditions. One biologically and physiologically relevant aspect of the HSR is the development of tolerance, i.e. the ability of a cell to become resistant to heat stress following a mild exposure to heat shock [70]. To understand this tolerance, Peper et al. [23] employed a basic model of HSR dynamics, which was in good agreement with observations in Reuber H35 rat hepatoma cells with both single and double heat-shock stimuli. Petre et al. [74] constructed a more detailed model incorporating components describing HSP-regulated transcription. They performed a sensitivity analysis to minimize the model by identifying and eliminating model components that only have marginal effects. Moreover, Scheff et al. [75] considered the effect of temperature changes on the transcriptional and multimerization processes. The above models could not capture data sets in different experimental conditions, yet after extension of the model by Petre et al. [74], a wide spectrum of data sets from HeLa cells were described well. From a local sensitivity analysis at steady state, this work also showed that the initial amount of HSF is expected to affect the stress-induced damage outcome much more than the amount of initial HSP and denatured proteins.

One of the HSFs, i.e. HSF2, is a key factor in developmental processes such as oogenesis, spermatogenesis, and corticogenesis [76]. HSF2 can cooperate with HSF1 by forming active heterotrimers in response to heat shock. A recent computational study employed probabilistic formulations revealing that the ratio of two HSF2 isoforms arising through alternative splicing controls the HSF1 transcriptional activity within individual cells [77]. However, no dynamical model has incorporated these recently proposed heterotrimeric processes, which might be explored when time-course data sets of the HSR become available during different HSF2-related developmental phases. As most models have been applied to experimental settings with thermal perturbations, linking of these models to drug-induced stress via HSP pathways needs to be further explored for its relevance to DILI in the future.

2.5. DNA damage response pathway

DNA damage can be caused both by endogenous processes, like deamination or replication errors, and by exogenous processes, like radiation and drugs. Drugs can cause DNA damage either directly or indirectly: Some drugs (mostly chemotherapeutics) directly intercalate between or covalently bind to DNA [78]. Other drugs lead to ROS that can bind to DNA and cause DNA damage, thereby indirectly activating the DNA damage response pathway [33,79–82]. When the DNA damage is not properly repaired, it can ultimately lead to cancer. Due to the relevance in cancer research, this pathway has been investigated in great detail, which includes modeling studies with the stressor γ-irradiation [25,83–97].

2.5.1. Mechanism

There are many different types of DNA damage, ranging from simple base substitutions to double-strand breaks (DBs). In this stress pathway, some form of DNA damage is the ‘controlled variable’ that activates the system, and activation of the pathway will lead to DNA repair. Because different types of damage need a different type of repair response or cell fate, the DNA damage response pathway consists of multiple modules, each dealing with different types of DNA damage [82]. One of the key elements (the transcription factor) of the pathway is the tumor suppressor p53 (Figure 2E), yet its activation can lead to different outcomes: The cell will try to repair the DNA damage, but to prevent further damage, it can also lead to permanent cell growth arrest (senescence) or apoptosis [98].

2.5.2. Modeling

Both p53 and mouse double minute 2 homolog (Mdm2), a negative regulator of p53, oscillate after getting stressed by DBs caused by γ-irradiation, both on population [83] and on single-cell level [99]. Lev Bar-or et al. [83] developed a minimal model of the interaction between p53 and Mdm2 to describe these oscillations. The transcription factor p53 stimulates the production of Mdm2, while Mdm2 enhances the degradation of p53, creating a negative feedback loop leading to an oscillation of both proteins when DNA damage is present. The oscillatory behavior only comes to a halt when the damage is repaired or is so excessive that the cell goes into apoptosis. Lev Bar-or et al. [83] argue that high levels of p53 are needed to activate the downstream DNA repair targets in adequate levels, yet continuously high levels of p53 induce apoptosis. Thus, oscillations give the cell the best chance to survive [83].

Not only p53 and Mdm2 oscillate after γ-irradiation, but also the active forms of the upstream regulators ataxia telangiectasia mutated (ATM) and checkpoint kinase 2 (Chk2) pulsate, and these pulses are required for p53 oscillations. Batchelor et al. [85] therefore proposed a model in which ATM and Chk2 activate p53 which subsequently activates (besides Mdm2) wild-type p53-induced phosphatase 1 (Wip1). Wip1 provides negative feedback in the pathway by inhibiting ATM and Chk2. The origin of the pulsation of both ATM and Chk2 is ongoing DNA damage, which keeps reactivating ATM and Chk2 after deactivation by p53 [85]. In a follow-up study, Batchelor et al. [84] showed that UV light instead of γ-irradiation also triggers p53 but is then regulated by different upstream targets: Whereas γ-irradiation triggers ATM and ATM triggers Chk2, UV light triggers ataxia telangiectasia and Rad3-related protein (ATR) and ATR triggers Chk1. Because ATR is not inhibited by Wip1, this results in a different response curve: after stimulation with UV light, p53 pulses only once [84].

Ma et al. [91] constructed a model that also included DNA repair itself i.e. consisting of three elements: a module about DNA damage and repair, a detector module, and an oscillator module. DNA damage and repair were modeled in a stochastic fashion, because the number of DBs at a low γ-irradiation level is very low. Specifically, DNA damage was described using a Poisson distribution based on the γ-irradiation level. Similarly, repair of DBs by DNA repair complexes was described by incorporating the amount of time needed to fix a DB as a stochastic variable. The DNA repair complexes also activate the second module, i.e. the detector of the system (ATM). ATM autoregulates itself in a positive feedback loop and activates the oscillators p53 and Mdm2. The added stochasticity ensures that there is some variability between individual cells. As a result, this modeling [91] showed that ongoing oscillations within individual cells can translate into the damped oscillations at the population level as observed experimentally. Moreover, variability in pulse number between cells has been confirmed in in vivo experiments [99].

An important cellular decision upon DNA damage is the fate of the stressed cell. When damage can still be repaired, the cell temporarily arrests the cell division process, and when the damage is too severe, the cell either goes into senescence or apoptosis. In cancer cells, p53 is often mutated, so despite substantial DNA damage, cells can continue to grow and divide. Some modeling efforts have been undertaken to integrate the DNA damage response dynamics with the decision to repair the DNA damage or to go into senescence/apopotosis. For example, Dolan et al. [25] successfully integrated two models, one that simulated senescence [92] and one that focused on damage repair [88]. The combined model captures a mechanistic overview of the DNA damage repair response, p53 signaling, and cell fate and predicts which dose of γ-irradiation leads to senescence. One of the clinically relevant findings was that repeated low doses of irradiation are as effective as a single high dose to steer cells into senescence, although the model overpredicted the effect of repeated low doses [25]. To our knowledge, p53 modeling has not been specifically applied to the liver and thus represents an important next step in DILI research.

2.6. Apoptosis pathway

Apoptotic cell death is a form of programmed cell death. It makes sure cells are destroyed in an organized manner, without losing the structure of the surrounding tissue [100] or bringing harmful molecules in its environment that could trigger inflammation [21].

2.6.1. Mechanism

Caspases are the primary agents of apoptotic cell death. They are a family of proteases that are inactive under homeostatic conditions. Caspases can be categorized in two groups, initiators and effectors. Initiator caspases are responsible for activating other caspases and amplifying the detected signal (i.e. the sensor of the system), while effectors are responsible for the dismantling of the cell. Initiation of the caspase pathway happens through either cell-extrinsic or -intrinsic stimuli. Cell-extrinsic signaling originates from outside the cell and is picked up by ‘death receptors’ located on the cell surface membrane, which dimerize and activate initiator caspase-8 or -10, which in turn activates effector caspase-3. The cell-intrinsic pathway starts in the mitochondria with mitochondrial outer membrane permeabilization: cytochrome c release from the mitochondrial intermembrane space into the cytoplasm leads to dimerization and activation of initiator caspase-9, which as before activates the effector caspase-3 [100]. The intrinsic apoptosis pathway can be activated by a wide variety of internal stressors, some discussed above like DNA damage and ER stress [101]. Therefore, the intrinsic pathway is the most important death pathway for DILI. Because apoptosis pathway modeling has already been often and extensively reviewed [101–108], we here focus on the highlights and its relevance to DILI.

2.6.2. Modeling

Several modeling studies have increased our understanding of apoptosis. Fussenegger et al. [109] published a model that included both cell-extrinsic and intrinsic death signaling. They used the model to test the potential for therapies inhibiting the apoptosis pathway, specifically applied to patients that have diseases in which cells are hypersensitive to apoptosis like Alzheimer (note though that caution should be taken when tempering with the apoptosis pathway because it safeguards us from for example cancer development). In the simulations, procaspase blocking was shown to be effective to prevent death through external signaling stimuli, but in silico, cells could still die from cytochrome c release due to internal signaling. Thus, because the cell has two distinct pathways, both should be blocked to guarantee that cells do not die from apoptosis. The model analysis further demonstrated that there exist a large number of redundant control mechanisms in the cell-intrinsic, but not in the cell-extrinsic caspase pathway. The general lack of control of the cell-extrinsic pathway may imply that this is in fact a ‘death override,’ which could explain why perfectly healthy cells go into apoptosis [109].

Eissing et al. [110] analyzed the bistability of the caspase cell-extrinsic pathway, because this is a mandatory property of the apoptotic pathway: On the one hand, cells should be resistant against minor fluctuations in signal in order not to trigger apoptosis by accident. On the other hand, when a cell is prompted into apoptosis, it should not be possible to revert back. The simplified model of the extrinsic pathway consisted of only six ODEs and considered stress to activate caspase-8, which subsequently activates caspase-3. Activated caspase-3 activates the remaining inactive caspase-8 molecules in a positive feedback loop. Moreover, active caspase-3 activates inhibitor of apoptosis (IAP, a member of a family of inhibitors of executioner caspases), which subsequently binds to active caspase-3 after which both get targeted for degradation. In their bistability analysis of this simple model, the ‘alive’ stage of the cell was unstable for biologically reasonable parameter values, implying that cells would go into apoptosis with only a little bit of noise. Eissing et al. [110] therefore suggest that caspase-8 needs a similar negative feedback loop as caspase-3 has via IAP. In the case of caspase-8, this role could be played by bifunctional apoptosis regulator (BAR), i.e. active caspase-8 activates BAR that binds to active caspase-8 after which both get targeted for degradation. These new model assumptions were sufficient to obtain robust bistability [110].

Chaves et al. [111] combined an NF-κB model [112] with the above-mentioned caspase model [110]. Because the connectivity between the NF-κB pathway and the caspase pathway is not completely understood, they added two known connections and tested various suggested but unconfirmed connections between the two models. The two know connections involve targets of NF-κB, i.e. NF-κB regulates both IAP and FLICE-inhibitory protein (FLIP). IAP as mentioned above inhibits active caspase-3, while FLIP inhibits an upstream target of caspase-8. When simulating the combined model with the two known links, it was apparent that a feedback mechanism was needed, to downregulate NF-κB upon caspase activation. Three feedback loops back from the caspase pathway to the NF-κB pathway were tested in the simulations. Inhibition from active caspase-3 on NF-κB indirectly via IKKα was most consistent with experimental observations. Thus, this analysis demonstrates how modeling can help to understand the links between stress (and death) pathways [111].

Linking activity of stress pathways to death pathways is an important step to obtain valuable insights into DILI. Cell death could be used as a risk indicator of the harmful effects drugs can have on the liver. Some of the above-discussed examples taking a first step in this direction (though not applied to the liver) include the prediction of cellular senescence upon DNA damage [25] and the effect of ER stress [24] or NF-κB activation [111] on apoptosis. Ultimately, different types of drugs and the associated activation patterns of diverse stress pathways need to be linked to cytotoxicity and to liver damage.

3. General discussion

Toxicity frequently occurs in the liver, because this organ is responsible for breaking down many xenobiotic compounds. Although many models have been developed for toxicity pathways, currently there are only a few models specifically designed for liver cells. As discussed, Leclerc et al. adapted a general oxidative stress model to describe liver cell dynamics [22]. The question arises whether all dynamical models designed for a particular organ or cell line can be applied to a different organ, in our case the liver. The underlying mechanisms of the stress pathways are likely to be highly similar amongst most cell types. The basic structure of the model can therefore be applied to different organs and cell types. However, some underlying system parameters, which are considered to be constant within a model, might differ between cell lines or cell types and so should be reevaluated. For example, Basak et al. [113] argue that their inflammatory stress pathway model can be adapted to other cell lines, after which some parameters need to be reset. Parameters based on the biophysical characteristics of the molecule such as the binding efficiency are unlikely to differ between cells, while other parameters can be affected by cell type-specific processes [113]. For example, the transcription rate of a gene could differ between cell lines due to chromatin regulation. As a second example, consider a protein (X) that strongly affects the binding of two other proteins (Y and Z), but whose expression is itself not modulated by the pathway. The pathway dynamics of two cell lines with either low or high expression of protein X can then be described by a model that does not explicitly consider protein X, albeit with different parameters for the interaction between proteins Y and Z. In conclusion, having a model that already describes the mechanism of the stress pathway of interest, although designed for very different cells, typically can be reused for other organs. Nevertheless, some parameters will need to be measured or fitted again. As already mentioned above, Zambrano et al. [48] showed that their simple model of NF-κB signaling could reproduce different types of oscillations, depending on values of their system parameters. They argue that the different behaviors observed in in vitro experiments are likely the effect of using different cell lines. This indeed is consistent with the hypothesis that adjusting the parameter values is often sufficient to reuse a model for a different tissue or cell line. The oxidative stress model by Leclerc et al. [22] focuses on the effect of a single drug and its mechanism of action. To be able to evaluate if drugs are potentially harmful to the liver, we need a model that is sufficiently general to be employed for many different drugs. Thus, such a model should take into account that drugs can activate a stress pathway in different ways. For example, the ER stress model by Trusina and Tang [65] was used to study the effect of different types of ER stress (translational and bond-breaking stress) on the dynamics of the UPR. Similarly, for the case of oxidative stress, metabolites can in some cases directly bind to cysteines on Keap1, or they can lead to increased ROS formation in mitochondria [27]. Moreover, different drugs preferentially bind to different cysteine sites on Keap1 and might therefore have a different effect. In conclusion, for the evaluation of drugs that are potentially harmful to the liver, we need models that are specifically adapted for the liver and that are applicable to different drugs.

To be of direct relevance to DILI, in the end, predictions by pathway models need to be made at an organ rather than cellular level. Insights gained by dynamical pathway models applied to the cellular level based on cell cultures may well be applicable to the human liver, but this is not necessarily straightforward because liver cell cultures are not fully representative of in vivo cellular behavior [114]: First, metabolism of drugs is typically poorly developed in cell lines compared to the in vivo situation, which may be problematic in case the metabolized form rather than the parent drug is the harmful compound. Second, multiple cell types may have a role in the development of adverse effects, e.g. during an inflammatory response. Solutions for these problems are being developed: Currently, in vitro toxicity tests are often done within (2D) cell monolayers. Newer methods apply toxicity tests on cell spheroids within 3D matrix gels and these have an improved drug metabolism compared to their 2D counterparts [114]. Progress is also made in cell cultures consisting of multiple cell types [115], and organ-on-a-chip approaches represent another promising solution [116].

One of the important issues when developing an ODE-model is the model size. Researchers often include a lot of detail of the signaling pathway they have studied. Basak et al. [113] have argued that the reasons for misregulation of a pathway and also potential therapeutic targets are often outside of the modeled part of the pathway. Therefore, in order to find therapeutic targets, larger models may indeed need to be studied. However, an overly complex model with many ODEs and parameters requires complicated analysis techniques and may lead to overfitting, whereas simpler models designed from first principles may just as well be able to describe experimental observations [48]. Moreover, for the specific application of DILI, small models may be sufficient to find a threshold at which potentially harmful compounds cause damage [113].

One option to deal with model complexity is to test new network hypotheses by studying subparts of a larger model. A good example of this approach is a study by Hofmann and coworkers on their NF-κB modeling [40]. Specifically, they tested different hypotheses on NF-κB homo- and heterodimer formation, initially with a very minimal model. They expanded the model until it explained the experimental observations on for instance dimer abundances. In the end, they altered their already well-established complex model, to include their new findings on dimer formation, and see if it still exhibits the same dynamics [40]. Another option to deal with model complexity is to first construct a complex model and employ reduction algorithms (e.g. [49]), which helps to determine the key processes involved in shaping the dynamics.

As discussed, cross-links between stress pathways and death pathways exist and are important to investigate in order to understand the relation between stress pathway activation and cell fate. One approach to obtain more insight is by linking models for these different pathways, but this also adds to the model complexity. Apart from cross-links between stress and death pathways, also the various stress pathways themselves are intertwined. For example, there is a suspected link between the oxidative and inflammatory stress pathways. One of the target genes of Nrf2 produces heme oxygenase 1 (HMOX1). HMOX1 upregulation can lead to inhibition of NF-κB-mediated transcription [117]. Indeed, activation of Nrf2 by some drugs (e.g. carbamazepine) inhibits the NF-κB response, yet not all drugs (e.g. acetaminophen) that activate Nrf2 have this effect on the NF-κB response [118]. The HMOX1 gene is also regulated by the transcriptional repressor BTB domain and CNC homolog 1 (BACH1) that inhibits under unstressed circumstances the binding of Nrf2 to the HMOX1 gene [119]. Thus, HMOX1 and BACH1 likely represent a selective cross-link between two of the discussed stress pathways.

Investigating the relations between stress pathways is needed for a full mechanistic understanding of adverse effects. However, a complicating factor is that cross links are not always activated and might be compound specific. We propose two methods to overcome this problem. As a first approach, the activation of all stress pathways upon exposure to a specific drug can be measured in vitro to determine which stress pathway models to use for this drug. This can for example be done in cell lines with built-in fluorescent reporters for the sensor, transcription factor, and downstream targets during live microscopy imaging [11]. For example, Wink et al. showed that known DILI drugs are frequently stress pathway specific, but some induce multiple stress pathways at a molecular level. The second, more challenging approach to take pathway cross talk into account is to create a general model that includes multiple relevant pathways and cross-links based on data in which the expression of proteins that are suspected to have a role in the cross talk is altered. However, this is an ambitious effort, which needs to start with the construction of suitable models for the individual pathways.

4. Expert opinion

The combination of in vitro testing of drugs on human cell lines, dynamical pathway modeling, and PK and PK–PD modeling have the potential to become a relevant screening method in order to detect DILI compounds at an early stage of drug development: In this procedure, cell lines are first used to assess the effects of compounds on relevant stress pathways within cells and cell fates such as death and proliferation arrest. This information is exploited to optimize and validate dynamical models describing these pathways. Subsequently, the models describe the effect of a compound on different pathways and cell fate in order to subsequently predict the tolerable dose. Finally, PK and PK–PD modeling are employed to translate this threshold to a corresponding tolerable drug intake in humans. This process was illustrated by Leclerc et al. [22] who integrated their dynamical model of the oxidative stress pathway with a PK model predicting how much of the drug actually reaches the cell. Experimental measurements were performed with an organ-on-a-chip setup, and the PK model described the amount of the drug flowing into the extracellular compartments, taking into account adsorption and desorption rates. Both the PK model and the oxidative stress model included a system variable to represent the amount of drug within the cell, which was used to bridge the two models. This approach facilitated the quantification of the effect of drugs added to the organ-on-a-chip on the stress induced to cells. In the future, such a strategy can be applied to other DILI-related drugs and different stress pathways.

This procedure is likely to work especially well for drugs that cause intrinsic DILI because they are expected to have a threshold with a similar range for different patients. Drugs that cause idiosyncratic DILI are hard to predict because of their low incidence rate. It could therefore be that the proposed combination of in vitro screening and computational modeling will not be successful in predicting idiosyncratic DILI until we have a better understanding of why these drugs are harmless to the majority of people but cause such extreme injuries in a few. An additional problem is that it is unlikely that in vitro cell lines react in the same way as cells within a patient, e.g. because additional perturbations are absent that naturally occur in patients, like cytokines that trigger immune responses.

Nevertheless, in order to improve our ability to detect idiosyncratic drugs, we need to obtain a quantitative understanding of the mechanism by which these drugs interact with the various stress and death pathways involved in liver injury. For example, drugs may alter the functioning of a stress or death pathway instead of merely activating it. Although this may be acceptable under some conditions, exposure to additional stress factors may lead to an unpredicted effect. Thus, focussing on the detection of altered pathway dynamics might give valuable clues to predict idiosyncratic drugs. Dynamical modeling is expected to have a pertinent role in this research.

To classify potential harmful drugs, it is important to determine which stress and death pathway is related to DILI. In vitro data can be employed to screen which stress and death pathways are activated by specific drugs. Furthermore, they are needed to parameterize and validate dynamical models. To this purpose, one of the most important factors includes how the behavior evolves over time. Thus, having measurements throughout the experiment, and not just end point measurements, is essential. Currently available methods include the use of reporter cell lines, in which the stress or cell death pathways can be monitored over time. In contrast to end-point assays, live confocal microscopy gives the opportunity to measure single-cell behavior over time [11].

In conclusion, to classify potentially harmful drugs, it is important to determine which stress and death pathways are related to liver injury. In vitro testing approaches can be employed to screen whether drugs activate or alter stress pathways. The development and application of in vitro tests and a quantitative understanding of the mechanisms behind stress and death pathways through dynamical pathway modeling together offer a promising prospect for DILI research in the years to come.

Article highlights

• Current techniques to predict DILI are insufficient, leading to health risks and financial losses.

• Dynamical modeling can give insights into mechanism, behaviour of a stress (death) pathway and its relation to adverse effects.

• In the last decade, enormous progress has been made in developing dynamical models for stress and (death) pathways.

• Combining high content in vitro experiments with dynamical modeling could lead to better techniques to detect adverse effects of drugs.

This box summarizes key points contained in the article.

Declaration of interest

This will include all details of conflict of interest, consultancies, honoraria, advisory board participation. The usual blurb will go after any declaration details The authors have no other relevant affiliations or financial involvement with any organization or entity with a financial interest in or financial conflict with the subject matter or materials discussed in the manuscript apart from those disclosed.

Acknowledgements

We would like to thank Suzanna Huppelschoten, Steven Hiemstra and Luc Bischoff for insightful discussions and Frederic Bois for his useful comments on a previous version of the manuscript.

Additional information

Funding

This work was supported by grants from the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO, InnoSysTox, number 40-42600-98-14016, to B van de Water and JB Beltman) and from the European Research Council (ERC, EUToxRisk21, number 681002, to to B van de Water and JB Beltman).