Modeling the dynamics of human liver failure post liver resection

Liver resection is an important clinical intervention to treat liver disease. Following liver resection, patients exhibit a wide range of outcomes including normal recovery, suppressed recovery, or liver failure, depending on the regenerative capacity of the remnant liver. The objective of this work is to study the distinct patient outcomes post hepatectomy and determine the processes that are accountable for liver failure. Our model based approach shows that cell death is one of the important processes but not the sole controlling process responsible for liver failure. Additionally, our simulations showed wide variation in the timescale of liver failure that is consistent with the clinically observed timescales of post hepatectomy liver failure scenarios. Liver failure can take place either instantaneously or after a certain delay. We analyzed a virtual patient cohort and concluded that remnant liver fraction is a key regulator of the timescale of liver failure, with higher remnant liver fraction leading to longer time delay prior to failure. Our results suggest that, for a given remnant liver fraction, modulating a combination of cell death controlling parameters and metabolic load may help shift the clinical outcome away from post hepatectomy liver failure towards normal recovery.


Introduction
Liver regeneration is a unique repair mechanism underlying physiological recovery following hepatic injury and enables surgical treatment as a viable clinical intervention into liver disease, small for size liver transplant and live donor liver transplantation. The process of liver regeneration takes place via hyperplasia and hypertrophy of hepatocytes as well as other liver resident cell types to reconstitute the liver morphology and function. In this process, the liver parenchymal cells, that is, differentiated post-mitotic hepatocytes, respond to signals from non-parenchymal cells that are stimulated by the hepatic injury and re-enter cell cycle, leading to multiple rounds of cell proliferation, thus compensating for the lost tissue mass [1][2][3].
Liver resection is widely performed for hepatocellular carcinoma (HCC), metastatic colorectal cancer and benign liver disease [4]. Liver response to resection is also relevant in live donors for liver transplantation, where the remnant liver in the donor needs to regenerate after removal of a portion of the liver for transplantation in the recipient. Advances in surgical techniques and expertise, careful patient selection and post-operative patient care have resulted in better clinical outcomes in the recent years [5]. However, Post hepatectomy liver failure (PHLF) is not uncommon due to a variety of risk factors, including patient-based factors such as age, weight, diabetes [6,7]. Preexisting liver disease such as cirrhosis, cholestasis, steatosis result in impaired liver regeneration [8].
Surgery related factors such as excessive blood loss, portal vein hypertension, ischemia-reperfusion injury and sepsis can result in liver failure. One of the most important factor that determines PHLF is how much liver can be resected which depends on the functional capacity of the remnant liver fraction. The percentage of liver failure post hepatectomy lies in the range of 0.7 to 9.1% [9]. To improve the outcome post-surgery, it is important to be able to predict the liver regeneration response based on patients' pre-and post-operative assessment. Mathematical modeling of liver regeneration can play a significant role in predicting the potential for PHLF prior to the surgical intervention, thus enabling consideration of alternative interventions to prevent or reduce chances of PHLF.
In this study, we started with our previously developed computational model of liver regeneration [10] and fine-tuned the parameters based on liver volumetric data from patients that underwent liver resection [11]. The main objective of this work is to model the dynamics of the liver failure scenario following a surgical resection. We focused on the model parameters that control the cell death process and generated a virtual patient cohort by sampling across the parameter space involving metabolic load and cell death sensitivity. We employed these virtual patients to analyze different modes of potential response; normal recovery, suppressed recovery and liver failure. Our simulations revealed wide variation in the timescale of liver failure consistent with the range of observed timescales in human PHLF scenarios. Simulations indicate that liver failure can either happen instantaneously or after a varying delay. We analyzed the distribution of this delay in a virtual patient cohort and found that the remnant liver fraction plays a significant role in controlling the time of delay in liver failure cases. Specifically, our simulations suggest that lower remnant liver fraction can lead to faster timescale of failure, depending on a balance of a subset of intrinsic parameters. We analyzed the effect of these controlling parameters corresponding to the metabolic load and cell death sensitivities on the overall cell death post resection. Our results suggest that specific combinations of these parameters can lead to a reasonable match to the clinically observed timescales of PHLF.

Mathematical Model
The mechanism of liver regeneration involves increase in portal pressure and production of cytokines from the Kupffer cells (resident macrophages in the liver). The cytokines activate the hepatocytes resulting in the priming and initiation of hypertrophy of the cells. Once the hepatocytes are primed they become responsive to the growth factor released from extracellular matrix via the action of matrix metalloproteinases (MMPs), as a consequence of which the parenchymal cells advance into the cell cycle and replicate leading to hyperplasia. The mathematical model (Cook et al. [10]) considered for this work incorporates the above-mentioned processes using eleven ordinary differential equations representing the hepatocytes in the three phases of the cell cycle quiescent (Q), priming (P) and replicating (R) states. The model comprises of seven equations for the molecular regulation, three equations for the cell state balances and the last equation accounts for relative cell mass growth, that is, hypertrophy. The system of equations is considered lumped since we assume the process of regeneration to be homogeneous in the remaining liver mass. The model equations are: where, Here, k 1 . . . k 7 are constants defined such that the molecular species are at steady state under the normal functioning of the liver. The initial conditions for solving the above-mentioned system of equations are: The system of equations for this model are stiff and were solved in Matlab using ode15s. The initial conditions for the quiescent state and molecular species signify the steady state under normal functioning of liver in which both primed and replicating cells are at zero levels.

Defining Virtual Patient Cohort
We employed the emerging approach of simulating a cohort of virtual patients [12,13] to assess the distribution of potential responses, in an attempt to account for the clinically observed wide variation of outcomes. In the present study, the virtual patients were defined based on sampling key parameters within a range of 2X, 3.125X and 7X, respectively, from the nominal values for the three controlling parameters (M, β cd , θ cd ). The sampling was performed in an unbiased approach to span the full range of the three-dimensional volume using Sobol sampling [14,15] to generate a cohort of 9000 virtual patients. The range of the controlling parameters were so chosen to capture the different response modes in the cohort of virtual patients under study.

Identification of Distinct Response Modes
The simulated results were classified into distinct response modes based on the level of liver mass fraction at the end of 2.5 years post liver resection-Normal recovery: 0.9-1.1, Suppressed recovery: <0.8. Patients are classified as exhibiting a liver failure scenario if the liver mass fraction dropped below 0.1 by 2.5 years post liver resection. The liver failure cases were further classified as exhibiting delayed failure or immediate failure based on the time taken to undergo failure. Immediate liver failure was deemed as occurring if the liver mass fraction is below 0.1 within the 5th day post resection and the remaining liver failure cases were classified as corresponding to the delayed failure scenarios.

Parameter Optimization
All the 33 parameters of the model (Table A1) were optimized using fmincon in Matlab using the elastic net approach with the regularization weightage of 0.001 for both Ridge and Lasso [16,17]. We explored the effect of using different regularization weights over a wide range from 0.0001 to 0.1. However, some of the regularization weights resulted in an integration error or terminated with a premature solution for a subset of the patients. Following this approach resulted in patient-specific regularization weights and thus likely yielded an overfitted model in some cases. We found that the regularization weight of 0.001 worked for all the patients and thus avoided over-fitting of the model.
The liver fraction at any given time point at and post resection is calculated as the ratio of remnant liver volume to the total liver volume prior to resection, based on the volumetric data from Yamamoto et al. [11]. However, note that the present ODE model considers fractional liver mass and not volume, as we incorporated the tissue growth dynamics in the model. For the purpose of the study, we considered the liver volume fraction information from the clinical data set to be equivalent to the liver mass fraction calculated in the simulations, under a reasonable assumption that the density of the tissue does not alter significantly after resection. This also permits us to compare the response profiles across individual patients that may have wide variation in the liver size. The model dynamics are sensitive to the three controlling parameters; metabolic load (M) and cell death sensitivities (β cd , θ cd ). The initial guess for the parameter optimization was considered from the parameter values of Cook et al. [10] for humans, except for the values of the three controlling parameters (M, θ cd , β cd ). These values were fixed from a representative case of a virtual patient undergoing delayed liver failure (Table A2).

Clinical Dataset
The human liver volumetry dataset used for the present work was obtained from the published literature [11], which contains information about the liver volume post resection in 196 patients. The dataset was analyzed for volume changes over time and a subset of 7 patients that showed liver failure were considered for further analysis and matching to model dynamics.

Model Repeatability and Reproducibility
The model equations and parameters for various simulations presented in the manuscript are available in the main text and Appendix, respectively. The Matlab code for virtual patient cohort simulation and response mode classification is available in the supplemental file S1. The Matlab code for parameter optimization is available as supplemental file S2. The network model in Equations

Modeling the Range of Response to Liver Resection in a Virtual Patient Cohort
We started with a computational model of liver regeneration, previously developed by our group [10]. The mathematical model of Cook et al. [10] considers hypertrophy and hyperplasia in liver following resection [18], by accounting for hepatocytes in different phases of the cell cycle. The molecular aspects of this network model of liver regeneration are largely based on the dynamics of cytokine and growth factor pathways stimulated by the liver resection injury. Immediately post resection, the hepatocytes undergo priming in response to the cytokines released by Kupffer cells, activating the Janus Kinase (JAK)-Signal Transducer and Activator of Transcription (STAT) signaling pathway in the hepatocytes. Primed hepatocytes become responsive to the growth factors linking the priming phase to the cell-cycle progression. In the present computational model, response to liver regeneration following partial hepatectomy is triggered by metabolic load per unit liver mass. The network model is shown in Figure 1A and the corresponding systems biology graphic notation (SBGN) [19] diagram for the liver regeneration model is provided as supplementary information ( Figure S1). This model can exhibit distinct modes of response ranging from complete recovery of liver mass to liver failure. As an illustrative example, Figure 1B shows the regeneration profile of two virtual patients that underwent 1/3rd liver resection but showed opposite responses, that is, recovery versus failure. Phase portrait [20,21] for the same virtual patients is shown in Figure 1C, depicting the evolution of hepatocytes in different phase of the cell cycle over time. Depending on the model parameters which corresponds to the two different virtual patients, for the same initial remnant liver mass, there exist two different stable steady states corresponding to recovery and failure. We simulated the dynamics of liver mass in response to resection in a virtual patient cohort generated by varying key model parameters (see Methods). We analyzed the distribution of responses in the virtual patient cohort and classified the virtual patients into multiple response categories based on the liver mass outcome at the end of 2.5 years following resection. We analyzed the responses of 9000 virtual patients to unravel the mechanisms that lead to different response modes. The virtual patient cohort was generated by varying the metabolic load and the cell death sensitivities via Sobol sampling, while the remaining 30 parameters (out of 33) were held at the same values as that of Cook et al. [10] for the case of human. We chose to vary metabolic load (M) and cell death sensitivities (β cd ,θ cd ) from the 33-dimensional parameter space since these 3 parameters control the cell death process as incorporated in the network model (see Equations (1)- (3) and (13)). In this formulation, liver failure is considered as likely to occur in two ways: (i) either the liver regeneration does not commence at the rate necessary to meet the increased functional demand post resection, or (ii) even as the regeneration process is initiated, the cell death rate is sufficiently high as to result in liver failure [22].
Our simulation-based analysis suggests that the virtual patients can be categorized into four distinct classes based on the response to the resection: normal recovery, suppressed recovery, delayed failure and immediate failure. The regeneration profiles of different classes of response for 10% and 33.3% are shown in Figure 2A,B,D,E. We note that the suppressed patients did not show fast recovery immediately after surgery. In Figure 2C we show a representative patient data demonstrating recovery post-surgery from Yamamoto et al. [11]. From Figure 2F we observe that a subset of the virtual patients exhibited instantaneous liver failure, whereas others showed a delayed liver failure), demonstrating a good match in the timescale of the liver failure based on volumetric data from patients in Yamamoto et al. [11]. Simulated regeneration profiles of liver failure response show a sudden drop in fraction of liver either immediately post resection or after a delay of few days. However, such a sudden drop was not readily apparent in the clinical data. It is unlikely that patients undergoing liver failure would be subjected to liver volumetric analysis in the clinic. Also, in reality, the liver mass does not become zero as is seen in the numerical simulations, since the patient's body might succumb to the non-functioning liver even before such condition of extremely low levels of liver mass is reached. However, the drop in the liver mass in simulations is instructive on the dynamics of the process, informing the potential trajectories taken by a range of patients, consistent with the wide range of time delays observed in patients exhibiting liver failure [11].
undergoing liver failure would be subjected to liver volumetric analysis in the clinic. Also, in reality, the liver mass does not become zero as is seen in the numerical simulations, since the patient's body might succumb to the non-functioning liver even before such condition of extremely low levels of liver mass is reached. However, the drop in the liver mass in simulations is instructive on the dynamics of the process, informing the potential trajectories taken by a range of patients, consistent with the wide range of time delays observed in patients exhibiting liver failure [11].

Figure 2.
Regeneration profiles of virtual patients depicting two types of liver recovery (Normal and suppressed) for different levels of resection (A) 10% resection (B) 33.3% resection (C) Volumetric liver data from Yamamoto et al. [1] for patients that exhibited liver recovery. Response modes of virtual patients that underwent delayed and immediate liver failure for (D) 10% resection (E) 33.3% resection (F) Volumetric liver data for patients showing liver failure from Yamamoto et al. [1].
We visualized the distribution of different classes of response modes-Normal, suppressed, delayed failure and immediate failure-as a function of the three controlling parameters (Figure 3). These classes were represented, albeit in different proportions, for varying levels of resection (10% and 33.3%). Figure 3A-F show the projection of classes of virtual patients for increasing levels of cell death sensitivity (θcd) at 10% and 33.3% resection. For low cell death sensitivity (θcd), we observe that liver failure occurs at high levels of metabolic load (M) and cell death sensitivity (βcd) and corresponds to only delayed liver failure ( Figure 3A,D). For higher level of resection, we observe that the normal recovery space is enveloped by suppressed recovery space. As the cell death sensitivity θcd increases, it becomes apparent that the parameter regions where normal recovery and liver failure occur are separated by a narrow region corresponding to the suppressed recovery ( Figure 3A-F). For either level of resection, low metabolic load and cell death sensitivities led to normal recovery. This appears to be governed primarily by low cell death ( Figure 4). By contrast, high values of these controlling parameters result in liver failure. Figure 3G,H show the pattern of the parameter region  33.3% resection (C) Volumetric liver data from Yamamoto et al. [11] for patients that exhibited liver recovery. Response modes of virtual patients that underwent delayed and immediate liver failure for (D) 10% resection (E) 33.3% resection (F) Volumetric liver data for patients showing liver failure from Yamamoto et al. [11].
We visualized the distribution of different classes of response modes-Normal, suppressed, delayed failure and immediate failure-as a function of the three controlling parameters (Figure 3). These classes were represented, albeit in different proportions, for varying levels of resection (10% and 33.3%). Figure 3A-F show the projection of classes of virtual patients for increasing levels of cell death sensitivity (θ cd ) at 10% and 33.3% resection. For low cell death sensitivity (θ cd ), we observe that liver failure occurs at high levels of metabolic load (M) and cell death sensitivity (β cd ) and corresponds to only delayed liver failure ( Figure 3A,D). For higher level of resection, we observe that the normal recovery space is enveloped by suppressed recovery space. As the cell death sensitivity θ cd increases, it becomes apparent that the parameter regions where normal recovery and liver failure occur are separated by a narrow region corresponding to the suppressed recovery ( Figure 3A-F). For either level of resection, low metabolic load and cell death sensitivities led to normal recovery. This appears to be governed primarily by low cell death (Figure 4). By contrast, high values of these controlling parameters result in liver failure. Figure 3G,H show the pattern of the parameter region corresponding to the virtual patients that exhibited delayed failure for remnant liver fraction of 0.9 and 0.667 (Supplementary information Figures S2 and S3, Video S1, S2). The proportion of cases that exhibited delayed liver failure varied from 8.73% to 40.22% of the total number of virtual patients for remnant liver fraction of 0.9 and 0.667, respectively (i.e., 10% and 33.3% resection). At the same time, the proportion of cases that exhibited suppressed recovery changed from 1.63% to 37.25% and the proportion of normal recovery cases reduced from 37.23% to 8.17%, as the remnant liver fraction decreased from 0.9 to 0.667 (i.e., increase in the level resection from 10% to 33%). This shift in the response of different virtual patients with increasing resection implies that the remnant liver mass significantly regulates the qualitative outcome of liver resection surgery. In the next section, we focus on the changes in the cell death process as a function of the three key controlling factors involving metabolic load and cell death sensitivities. exhibited delayed liver failure varied from 8.73% to 40.22% of the total number of virtual patients for remnant liver fraction of 0.9 and 0.667, respectively (i.e., 10% and 33.3% resection). At the same time, the proportion of cases that exhibited suppressed recovery changed from 1.63% to 37.25% and the proportion of normal recovery cases reduced from 37.23% to 8.17%, as the remnant liver fraction decreased from 0.9 to 0.667 (i.e., increase in the level resection from 10% to 33%). This shift in the response of different virtual patients with increasing resection implies that the remnant liver mass significantly regulates the qualitative outcome of liver resection surgery. In the next section, we focus on the changes in the cell death process as a function of the three key controlling factors involving metabolic load and cell death sensitivities.

Cell Death as a Function of the Controlling Parameters in Different Classes of Patients
Following partial hepatectomy, both positive and negative stimulus are triggered resulting in liver regeneration or failure based on a balance of these stimulus. Cell death is likely an important regulatory process in liver regeneration, levels of which determine whether the organ can recover following injury. Analysis of the cellular equations (Equations (1)-(3)) of the model along with the cell death sigmoidal function σ cd (Equation (13)) reveals that there are three controlling parameters: metabolic demand (M) and the two cell death sensitivity parameters (θ cd ) and (β cd ) that influence decay of cells post resection. In this section, we analyze the cell death profile as a function (Equation (13)) of these three parameters (M, θ cd , β cd ), individually and in combination. A change in the cell death sensitivity β cd results in a shift of the angle or slope of the cell death function (Figure 4A). At a high level of resection (i.e., low fraction of remnant liver), increasing β cd results in a downward shift in the cell death function, lowering the cell death ( Figure 4A). However, at low level of resection (i.e., high remnant liver mass fraction) the cell death increases with increasing β cd ( Figure 4A). The cell death profile as a function of increasing cell death sensitivity (θ cd ) or metabolic load (M) results in a horizontal shift towards higher levels of remnant liver fraction ( Figure 4B,C). Consequently, for a given level of remnant liver fraction, increasing metabolic load (M) or cell death sensitivity (θ cd ) lead to higher levels of cell death. The cell death profile when both the cell death sensitivity (β cd ) and metabolic load are changed simultaneously show a wide variation in both slope and the switching threshold of the sigmoidal cell death function ( Figure 4D).

Cell Death as a Function of the Controlling Parameters in Different Classes of Patients
Following partial hepatectomy, both positive and negative stimulus are triggered resulting in liver regeneration or failure based on a balance of these stimulus. Cell death is likely an important regulatory process in liver regeneration, levels of which determine whether the organ can recover following injury. Analysis of the cellular equations (Equations (1)-(3)) of the model along with the cell death sigmoidal function σcd (Equation (13)) reveals that there are three controlling parameters: metabolic demand (M) and the two cell death sensitivity parameters (θcd) and (βcd) that influence decay of cells post resection. In this section, we analyze the cell death profile as a function (Equation (13)) of these three parameters (M, θcd, βcd), individually and in combination. A change in the cell death sensitivity βcd results in a shift of the angle or slope of the cell death function (Figure 4A). At a high level of resection (i.e., low fraction of remnant liver), increasing βcd results in a downward shift in the cell death function, lowering the cell death ( Figure 4A). However, at low level of resection (i.e., high remnant liver mass fraction) the cell death increases with increasing βcd ( Figure 4A). The cell death profile as a function of increasing cell death sensitivity (θcd) or metabolic load (M) results in a horizontal shift towards higher levels of remnant liver fraction ( Figure 4B,C). Consequently, for a given level of remnant liver fraction, increasing metabolic load (M) or cell death sensitivity (θcd) lead to higher levels of cell death. The cell death profile when both the cell death sensitivity (βcd) and metabolic load are changed simultaneously show a wide variation in both slope and the switching threshold of the sigmoidal cell death function ( Figure 4D).  We analyzed the cell death function of the 9000 virtual patients that were categorized for response modes in Section 3.1. The cell death functions corresponding to the four response mode classes of patients are shown in Figure 5. The cell death functions of both recovery classes, normal and suppressed growth, are steep and shifted towards lower remnant liver fraction ( Figure 5A,B,E,F), which resulted in almost negligible cell death at low levels of resection (i.e., high remnant liver fraction). The cell death functions for the liver failure classes show wide variation and are more inclined as compared to normal and suppressed recovery classes, resulting in higher cell death for any given level of resection (Figure 5C,D,G,H). The ranges of cell death functions are overlapping across the response modes, suggesting that cell death function is not the sole determining factor for classifying the fate of the patient post resection surgery. Specifically, the level of resection also plays a crucial role in determining the dynamics of liver recovery and potential failure.

Comparison of Cell Death between Recovery versus Failure Scenarios for Varying Level of Resection
In this section, we compare the responses of two virtual patients showing liver recovery and liver failure to different levels of resection to examine the distribution of responses. Figure 6A,D shows the regeneration profiles of a representative normal recovery and a liver failure cases from the virtual patient cohort for 1/3rd resection. The response profiles were simulated for different levels of resection. The virtual patient shows recovery for varying levels of resection, until a threshold value of 83%, beyond which liver failure occurred due to lack of a robust regenerative response ( Figure 6A). By contrast, the response of the patient with controlling parameters leading to liver failure was insensitive to the level of remnant liver mass ( Figure 6D). Simulations suggest that at higher levels of resection (i.e., lower remnant liver fraction), the patient can exhibit an inverse response in which the liver mass enters an initial recovery phase while undergoing a continuous loss of liver tissue that overcomes the recovery, resulting in liver failure. The failure is reflected in the response profile as a fast decline in the remnant liver fraction and the time delay of failure was longer with higher levels of remnant liver mass ( Figure 6D).
For the above two virtual patient cases, we examined the evolution of hepatocytes in quiescent versus proliferating phase over time for different levels of resection. The corresponding projections of the phase spaces shown in Figure 6B,E suggest the existence of two attractors, one corresponding to failure and the other to recovery. For the first virtual patient case, a threshold of 83% resection separates the two attractors, such that increasing level of resection results in a shift in the stability of the system from one attractor to another, that is, exhibiting a transition from recovery to failure ( Figure 6B). By contrast, the threshold of failure for the second virtual patient was low, likely corresponding to the existence of only one attractor for this case that is associated to liver failure ( Figure 6E). We compared the cell death function for the two cases for uncovering potential differences over time. The surface plots shown in Figure 6C,F correspond to the two virtual patient cases and the trajectories on the surface depict the cell death process for a specific level of resection. For the first virtual patient case, the cell death is negligible for low levels of resection as seen by the trajectories largely spanning the lower part of the cell death function (Figure 6C). At a resection level beyond the threshold of failure, there is a drastic change in the cell death profile with the trajectories distributed in the upper portion of the cell death function, resulting in high cell death. In the second virtual patient case, the cell death trajectories show initial fluctuations at low levels of the cell death function, before shifting upward and saturating to the maximum cell death which leads to a fast decline in the remnant liver fraction ( Figure 6F).

Tuning the Model for Different Liver Failure Patients
The analysis so far has provided insights into the effect of key biochemical and biophysical parameters in the liver regeneration process. We used this knowledge to tune the model parameters to capture the timescale of response to resection as reflected in the volumetric liver data for the patients that exhibited liver failure available in Yamamoto et al. [11]. Figure 7 shows the regeneration profiles predicted and optimized parameters for 7 real patients who exhibited liver failure. The model parameters were optimized for each patient-specific data set using elastic net (see Methods). The initial values for parameter optimization were chosen such that the values of the controlling parameters (M, θ cd , β cd ) were set from a representative case of a virtual patient belonging to the class of delayed liver failure and the remaining 30 parameters were set from Cook et al. [10]. The simulated liver fraction profiles from models with patient-specific optimized parameters showed a good match with the volumetric data and captured both the slow (ID114, 123, 131, 29) and fast (ID45, 245, 135) timescales of liver failure ( Figure 7A). Comparing the optimized parameters to the initial values show that the priming of hepatocytes and JAK-STAT pathway contribute to the difference in the timescale of liver failure ( Figure 7B). The table of optimized parameter values for the individual patients is given in Appendix A.
The profiles of remnant liver fraction for the slow and fast timescale of liver failure cases are shown in Figure 8A,F. The corresponding phase space portrait projection and the cell death functions of patients undergoing liver failure at the slow and fast timescales are shown in Figure 8A-E and Figure 8F-I, respectively. Even as the regeneration profiles of the patients with slow timescale of response appear to match a suppressed recovery case, the phase plane depicts that the patients likely undergo liver failure, as the fraction of remnant liver in each case is progressing towards the attractor of zero levels. Analysis of the patient-specific cell death functions and trajectories suggest that the cell death trajectories of patients showing slow timescale of failure were largely restricted to the lower levels of cell death function. By contrast, the cell death trajectories of patients with relatively faster timescale of liver failure continually progressed from lower levels with a steep ascent to reach maximal levels of cell death function, corresponding to a rapid decline in the predicted liver fraction.
We note that the optimized values of the controlling parameters of all the patients were similar (Appendix A) and yet the cell death trajectories were rather divergent across the patients (Figure 8B-E,G-I). This is because cell death is not only dependent on the identified controlling parameters but also on the remnant liver fraction, which depends on the regenerative processes including priming pathways and proliferation. The parameters for the regenerative capacity are different between the patients ( Figure 7B) and these differences manifest in the net amount of recovery, or lack thereof, in the fractional liver, indirectly affecting the cell death trajectory. Our simulations suggest that patients with slow timescale of liver failure sustain low cell death for a longer duration compared to the patients with a relatively faster timescale of failure that experience high cell death with a rapid decline in the remnant liver fraction leading to failure and patient death. These results from analysis of patient-specific model parameterization and simulation provide insights into how the cell death trajectories likely control the timescale of response to resection, particularly for the liver failure cases.

Discussion
The focus of this study is to understand the interplay between the different processes that are activated after a liver resection and can lead to post hepatectomy liver failure. We started with a computational model of liver regeneration process based on the published cellular network model of Cook et al. [10]. We simulated the model over a wide range of parameters and found that the model captures the wide variability in patient outcomes of liver regeneration and accounts for the variation in the observed timescales of liver failure. The variability in response to liver resection may arise due to inter-individual differences in the sensitivity to injury, which is likely due to disease etiology as well as intrinsic genetic factors [23]. For instance, hepatic steatosis alters the potential for cell survival and proliferation after liver resection as compared to the normal liver [10,24].
We identified the controlling parameters in the model which govern the timescale of the post hepatectomy liver failure in humans. These controlling parameters correspond to the metabolic load and cell death sensitivities. We generated a cohort of virtual patients by varying the controlling parameters and classified the virtual patients based on their model-predicted outcomes post liver resection. Our analysis suggests an inverse relationship between the levels of resection on the proportion of cases showing normal growth, suppressed recovery, with concomitant increase in the liver failure cases. These results are consistent with expectations of increased potential for failure with increasing level of resection in humans and animals [6]. Our model-based virtual patient analysis revealed that the parameter space corresponding to the transition from recovery to failure outcomes corresponds to delayed liver failure and the proportion of the delayed liver failure is dependent on the level of resection. These results show that the remnant liver fraction is an important decision making variable that not only discriminates between liver failure and recovery but also governs the timescale of failure especially when a patient experiences delayed liver failure. These results provide new insights in interpreting the wide range of variation seen in clinical data sets [11].
Our model-based results on the critical role of a combination of cell death process and remnant liver fraction in controlling the liver regenerative response is consistent with the findings from animal studies showing that beyond a certain threshold of resection, cell death increases substantially, resulting in high likelihood of failure [6,25]. Importantly, reducing the cell death rescues the animals from undergoing failure after partial hepatectomy [26]. Similarly, preventing the hepatocyte death by blocking mitochondrial permeability transition can improve liver regeneration even for low remnant liver fraction in animals [27,28]. Our model-based results place these experimental findings in the context of cellular networks that underlie the regeneration process, as well as point out additional parameters that can be potentially manipulated to help shift the trajectory of the response modes to the recovery zone. The key controlling parameters were determined such that their combination decided the level of cell death that occurred post resection. With our analysis, we concluded that though cell death is an important process which influences the liver regeneration it is not the sole factor which determines liver failure. Liver failure is conditional upon the load and other perioperative conditions. This is likely to be patient specific and needs to be correlated with epidemiological, clinical measures and disease etiology [29]. Modulation of these additional factors likely improve the probability of liver recovery and hence have the potential to yield additional clinical options to reduce the chances of liver failure due to lack of regenerative response.
Various parameters of this model were derived based on biochemical data and animal/clinical observations (Furchtgott et al. [30]; Cook et al. [10]). The validation relevant to our present study is regarding the dynamics of liver failure. Model simulations predict a wide range of time to failure, which is consistent with the variability of the time scale of liver failure, seen in human data (Yamamoto et al. [11]). We note that the virtual patient analysis approach differs from patient-specific modeling that identifies parameters for individual patients and then performs comparative analysis across patients. By contrast, the virtual patient simulations are evaluated against the types of behavior (e.g., recovery versus failure), distribution of outcomes (e.g., range of time scales of failure) and so forth., that are exhibited in a cohort (An et al. [12,13]). In the present study, information from the virtual patient analysis on delayed liver failure response was used for optimization of the model parameters for patient-specific data, which were not used to build the original model. Our analysis identifies cell death relevant parameters and metabolic load as key factors controlling the time scale of liver failure post-surgical resection. These results are consistent with Yamamoto et al. [11] findings from analysis of patient data that suggest the extent of blood loss during surgery as a key factor discriminating between failure versus recovery scenarios. It is likely that perioperative factors such as blood loss manifest their effects on liver at least partly through increased cell death and higher load on the remaining functional organ.
Our study was focused on the regulatory networks at the cellular and pathway scales in the liver during the response to resection. However, failure or success of the resection surgery also depends on the recovery of the metabolic function in addition to the liver mass. The present computational framework can be extended to incorporate metabolic components to account for relationship between functional liver fraction and key metabolic functions, for example, glucose homeostasis [31], ammonia metabolism and urea cycle [32], or genome-scale metabolic models [33]. Liver failure can also occur due to extra-hepatic factors including immune response, as well as inter-organ metabolic and physiological relationships. For example, hepatic encephalopathy is an important condition that is guarded against by monitoring ammonia homeostasis, which is considered as a clinically-relevant indicator of liver metabolic function after liver surgery [34]. These processes are not explicitly accounted for in the present network model. One can consider that metabolic load parameter serves as a lumped, phenomenological factor that accounts for the stress on liver, with metabolic load per unit of liver mass (M/N in Equations (4), (9) and (11)) serving as a stimulus for regeneration as well as cell death [10,30,35]. Opportunities exist for integrating the liver regeneration model with multi-organ physiology models to expand the utility of the model into a whole body context. The model-based research discussed in this work on liver regeneration has been based on a lumped model considering the hepatocyte functional states. This approach can be extended to incorporate functional states of other liver cell types including hepatic stellate cells, sinusoidal endothelial cells and Kupffer cells (resident macrophages of the liver) [36], as well as potential emerging rescue approaches such as stem cell transplant [37]. Taken together, the above detailed extensions to the network model will permit relating the key controlling parameters to inter-cellular, tissue-scale and whole body physiological scale parameters and will likely provide additional insights into novel venues for clinical management and intervention.
Our use of a virtual patient cohort approach helped delineate the key parameter intervals and combinatorial dependencies that correspond to a wide range of predicted outcomes post resection. Fruitful next steps could be to correlate these model-predicted parameter subspaces to patient demographical data as well as preoperative clinical information and disease etiology. Development of such a correspondence between patient information and model parameters is likely to aid in generalized application of the dynamic modeling to a wide range of liver surgery scenarios. Such a model-based approach informed by patient data to constrain the parameters can assist in clinical decision making by predicting the categorical outcome prior to the clinical intervention. Specifically, our model-based approach can aid in estimating the safe level of resection a patient can undergo along with predicting the need for necessary manipulations of cell death sensitivities and metabolic load to maximize the chances of recovery.