Computational Hypothesis: How Intra-Hepatic Functional Heterogeneity May Influence the Cascading Progression of Free Fatty Acid-Induced Non-Alcoholic Fatty Liver Disease (NAFLD)

Non-Alcoholic Fatty Liver Disease (NAFLD) is the most common type of chronic liver disease in developed nations, affecting around 25% of the population. Elucidating the factors causing NAFLD in individual patients to progress in different rates and to different degrees of severity, is a matter of active medical research. Here, we aim to provide evidence that the intra-hepatic heterogeneity of rheological, metabolic and tissue-regenerating capacities plays a central role in disease progression. We developed a generic mathematical model that constitutes the liver as ensemble of small liver units differing in their capacities to metabolize potentially cytotoxic free fatty acids (FFAs) and to repair FFA-induced cell damage. Transition from simple steatosis to more severe forms of NAFLD is described as self-amplifying process of cascading liver failure, which, to stop, depends essentially on the distribution of functional capacities across the liver. Model simulations provided the following insights: (1) A persistently high plasma level of FFAs is sufficient to drive the liver through different stages of NAFLD; (2) Presence of NAFLD amplifies the deleterious impact of additional tissue-damaging hits; and (3) Coexistence of non-steatotic and highly steatotic regions is indicative for the later occurrence of severe NAFLD stages.


Introduction
Non-Alcoholic Fatty Liver Disease (NAFLD) is the most common type of chronic liver disease in developed nations, where it affects around 25% of the population [1]. NAFLD starts with simple steatosis and may progress to severe liver diseases like cirrhosis or hepatocellular carcinoma (HCC). Up to 30% of patients with liver steatosis develop a Non-Alcoholic Steatohepatitis (NASH), an inflammatory state of the liver that in about 20% of cases develops further to cirrhosis and end-stage liver failure [1]. Owing to the high metabolic reserve capacity of the liver, progression of NAFLD typically proceeds silently over long periods before an abrupt worsening of central metabolic liver functions results in severe clinical symptoms requiring urgent treatment.
Steatosis of the liver in the absence of significant alcohol intake is commonly associated with insulin resistance resulting from hyper-caloric diet and a sedentary lifestyle. Insulin resistance of the adipose tissue increases the release of free fatty acids (FFAs) into the blood [2] and thus promotes the uptake of FFAs into liver, kidney, heart and other organs using FFAs as preferred energy-delivering substrates. FFAs exceeding the cellular need for ATP synthesis are esterified to triacylglycerol (TAG) and membrane lipids. Conversion of excess fatty acids into complex lipids serves as a detoxification mechanism [3][4][5][6] as elevated levels of FFAs and some of their reaction products may cause cell damage [7]. As shown in mice with downregulated diacyglycerol acyl transferase 2 (DGAT2), a rate-limiting enzyme of TAG synthesis, inhibition of TAG synthesis improved hepatic steatosis but exacerbated liver damage exerted by elevated levels of potentially cytotoxic FFAs [6]. FFAinduced hepatotoxicity includes various molecular mechanisms as, for example, induction of endoplasmic reticulum stress, generation of free radicals and subsequent activation of the mitochondrial apoptosis pathway [8]. Necrosis and apoptosis of hepatocytes mount a dynamic multicellular response wherein stromal cells are activated in situ and/or recruited from the bloodstream, the extracellular matrix is remodeled, and epithelial cells expand to replenish their lost numbers [9]. If the FFA challenge persists, a quasi-stationary balance between damage and tissue regeneration is established, which over the time may slowly shift towards fewer and fewer vital hepatocytes and more and more non-functional fibrotic lesions [9].
Several factors may explain why patients with diagnosed steatosis may or may not develop NASH and even more severe chronic liver disease. As with all diseases, genetic factors may determine the susceptibility of an organ to damaging events. Genome-wide association studies have revealed several single nucleotide polymorphisms associated with the pathology of NAFLD, among them the gene variant I148M of the enzyme Patatin-like phospholipase domain-containing 3 plays an important role in the cellular lipid and lipid droplet metabolism [10]. On top of this, the liver is continuously confronted with all kinds of orally administered toxins and gut-derived pathogens.
A now widely accepted hypothesis postulates that simple steatosis ("first hit") has to be followed by further pathophysiological hits, such as pathogen-associated acute inflammation or gut-derived endotoxins, to push the liver successively into a critical state [11]. A third factor influencing the progression of NAFLD consists in the intra-hepatic spatial heterogeneity of metabolic and immune-modulatory capacities. Gradients along the porto-central axis of the liver lobule exist not only for genes and proteins of hepatocytes leading to zone dependent differences in lipid metabolism [12] but also in other types of liver cells and the matrix of the space of Disse [13,14]. Zone-dependent differences in gene expression within a single lobule can be partially accounted for by concentrations gradients of metabolites, hormones and morphogens along the sinusoidal blood stream [15,16]. In adults, steatosis is most intense around the central veins (predominantly in zones 2 and 3), whereas children may have an alternate pattern of progressive NAFLD characterized by a zone 1 distribution of steatosis, inflammation and fibrosis [17].
In computational liver research, mathematical models can reveal hidden principles of organ function by classification of individual components, analyzing their interactions and simulating the effects that are difficult to investigate experimentally [18]. A recent modeling approach focused on molecular mechanisms accounting for differences in the steatosis susceptibility of portal and central regions of the acinus [19]. Another systems biology approach coupled in vitro experimentation to a multi-scale mathematical model to investigate the roles of various sugars and fat in NAFLD pathogenesis [20]. In this work, we use a multi-scale model to shed light on a novel aspect that hitherto has not been implicated in NAFLD progression: the functional heterogeneity of the liver on a spatial scale that is far beyond the intra-lobular dimension. We refer to this type of heterogeneity as "macroscale heterogeneity" as it runs on lengths scales of millimeters and centimeters, which are orders of magnitude above the size of a single liver lobule. Support for the existence of macro-scale heterogeneity comes mainly from functional liver imaging and histological findings. High-resolution imaging of hepatic blood flow revealed large regional variability up to a factor of three [21,22]. Hepatic clearance of 2-[ 18 F]fluoro-2-deoxy-D-galactose measured by PET/CT varied by about 24% in patients and 14% in healthy subject [23]. The spatial distribution of liver fat in adults with NAFLD showed a variability in the range of 0.7%-4.5% [24]. Patient-specific spatial pattern of fat dispositions may largely vary from diffuse fat accumulation with and without focal sparing and focal fat accumulation [25]. Presence of macro-scale heterogeneity in functional and metabolic parameters is also reflected by differences of NAFLD-associated histopathologic lesions between different anatomical parts of the liver [26].
Taken together, progression of NAFLD can be affected by several intrinsic and environmental factors, among which the rise of the serum FFA level, the distribution of intrahepatic functional capacities, and the eventual occurrence of additional liver-damaging events play a central role. The question of how these factors may act together in the patient-specific development of NAFLD can hardly be addressed in a clinical approach for obvious technical and ethical reasons. Computer simulations based on physiologically reliable mathematical models may serve as an alternative in such situation. This prompted us to develop a generic and in many ways further expandable mathematical model of NAFLD progression.

Materials and Methods-Mathematical Model
The mathematical model used for our disease simulations consists of three modules: (1) The hemodynamic module describes the blood flow within the liver sinusoids as was previously described [27]; (2) the metabolic module describes the cellular turnover of FFAs and TAG; and (3) the damage-repair module describes the FFA-induced damage of hepatocytes and tissue repair.

Module 2: Kinetic Model of FFA and TAG Turnover in a Liver Unit
We subdivided the liver into a larger number of small liver units (LUs) differing in their hemodynamic properties, metabolic capacities, vulnerability against FFA-induced cell damage and wound-healing capacities. A single LU is described by a two-compartment model: the vascular compartment formed by the sinusoids receives fatty acids with the blood from the hepatic arteries and exchanges fatty acids with the cellular compartment, which is equipped with the average metabolic capacities of all hepatocytes in the unit (see Figure 1). The numbers on the reaction arrows in Figure 1B are FFA particle numbers (in micromoles) and particle transport rates (in micromoles/100 mL/min) for a healthy LU with a volume of 100 mL. They represent typical values of transport and metabolization rates of FFAs measured in the healthy human liver (see references in Table 1). in particular the skeletal muscle (utilization rate v ET ). (B) FFAs enter the vascular bed of the liver with the rate v ps that is determined by the plasma concentration of fatty acids and the blood flow rate. A certain fraction of fatty acids is taken up by hepatocytes with rate v u , the remaining part re-enters the circulation via the venous blood efflux. Cellular fatty acids (FFA c ) are either formed de novo (synthesis rate v denovo ) or taken up from the plasma. They can be used as building blocks for the synthesis of various lipids, mainly triacylglycerols (TAG) (synthesis rate v +TAG ). Note that the variable TAG denotes the pool of FFAs esterified in TAG in the model. Alternatively, FFA c can be degraded to acetyl-CoA, either by mitochondrial and peroxisomal β-oxidation or CYP-450-mediated ω-oxidation (oxidation rate v ßox ). TAG, the most abundant cellular lipid, is stored in lipid droplets, which can by hydrolyzed to FFA c if needed (degradation rate v -TAG ). A fraction of TAG is loaded in VLDL lipoproteins and exported (export rate v VLDL ). The numbers on the reaction arrows are FFA particle numbers (in micromoles) and particle transport rates (in micromoles/100 mL/min) for a healthy LU with a volume of 100 mL. The corresponding mass transport rates are given in Table 1. Table 1. Rate laws for lipid fluxes and values of kinetic parameters for a standard liver unit. Rate constants were based on the particle numbers and particle transport rates shown in Figure 1. Numerical values of the rate constants were chosen such that for a healthy liver under steady-state conditions, the metabolic rates of the model were in good concordance with reported average experimental data (column 5). In the rate laws v ap and v vp for the arterial and venous plasma flow, the parameters r s (3.15 µm), l s (375 µm) and Ω s (1.169·10 −11 L) denote the average radius, length and volume of sinusoids in the LU. For details, see Berndt et al., 2018 [27]. *) calculated from the published VLDL-TAG secretion rate [28] by assuming a liver volume of 1350 mL.

Flux
Meaning Blood flow rate of 1 mL blood/1 mL liver The kinetics of FFAs in the plasma (FFA p ), in the sinusoid (FFA s ) and in the hepatocyte (FFA) in the i-th LU (i = 1, . . . ,N LU ) is governed by the following system of first-order differential equations: The first equation describes the concentration change of plasma FFAs as the net outcome of three processes: Release of FFAs from the adipose tissue with rate v FFA at , the utilization of FFAs by extra-hepatic tissue (predominantly skeletal muscle and kidney) with rate v FFA et , and the net uptake of FFAs by the liver given by the arterio-venous concentration difference multiplied by a scaling factor, which relates the blood volume received by the liver to the total blood volume. The second equation describes the concentration change of FFAs in the sinusoidal blood stream. The parameter η i , 0 ≤ η i ≤ 1, denotes the fraction of metabolically intact hepatocytes in the i-th LU. The third equation describes the concentration change of FFAs in the hepatocytes as the result of five elementary processes: Uptake from the external space, mitochondrial β-oxidation to acetyl-CoA, esterification to TAG, hydrolysis of TAG to FFAs and glycerol, and de novo lipid synthesis. The quadratic term in the rate equation v +TAG for TAG synthesis takes into account that with increasing cellular content of FFAs, the synthesis of TAG and storage in lipid droplets proceeds super-linear due to the upregulation of lipogenic enzymes [37]. FFA crit denotes the cellular concentration threshold of fatty acids at which the rate of TAG synthesis starts to rise in a non-linear fashion. The fourth equation describes the concentration change of cellular TAG due to the formation and hydrolysis of TAG and the packing into VLDL lipoproteins and their export to extra-hepatic tissues. The fifth equation describes the removal of fatty acids from the plasma by extrahepatic tissue.
The blood flow is given by: Again, the upper index (i) numbers the LUs and r s i , l s i , Ω s i denote the radius, length and volume of the i-th LU. vis blood is the blood viscosity, ∆p a i is the pressure gradient between arterial blood and the i-th LU and ∆p v i denotes the pressure gradient between venous blood and the i-th LU. For details, see Berndt et al., 2018 [27]. Figure 2 illustrates how differences in the metabolic capacities of individual LUs may influence the diurnal level of cellular FFAs. The rate of fatty acid influx into the plasma was modeled by a periodic function with period T, reflecting the insulin-dependent variation of FFA release from the adipose tissue to the plasma.  Table 1. (B,C) Impact of the cellular TAG storage capacity on variations of cellular FFAs and TAG. The curves were generated by varying the rate constant k +TAG for TAG synthesis in small steps between 0.036 min −1 and 3.6 min −1 . The values of all other model parameters were the same as given in Table 1. Red curves: k +TAG = 0 (no TAG synthesis). Green curves: k +TAG = 3.6 min −1 (highest TAG synthesis). (D,E) Time-course of FFAs and TAG at randomly varied kinetic parameters. One hundred curves were generated with random, normally distributed parameter values having as sample mean the original values in Table 1 and possessing a standard deviation of 20% (=0.2 mean value).
The amplitude of variations of the lipid species decreased in the order FFA p → FFA s → FFA cell → TAG. Figure 2B,C illustrate how daily concentration changes of cellular FFAs depend on the TAG-synthesizing capacity. Increasing the cellular capacity to convert FFAs into TAG reduces the amplitude of diurnal variations and the 24 h mean of the cellular FFA concentration, but increases the amplitude of variations and the mean of the TAG content. Note that at normal rate of TAG synthesis the daily maximum of cellular TAG is about 1.8-fold higher than the minimum (cf. green curve in Figure 2C). This difference is in good agreement with measured TAG variations between the fasted and the fed state [39]. Figure 2D,E illustrate the effect of random variations of metabolic parameters on the cellular profile of FFAs and TAG. Twenty percent of random variations of the reference parameters (given in Table 1) may give rise to concentration maxima of cellular FFAs between 0.08 mM and 0.68 mM. For the variation of TAG, the coefficient of variation (CV; standard deviation to mean) was 35%, which is in good concordance with values obtained in outbred liver-healthy Wistar rats [40].

Module 3: Fatty-Acid Induced Cell Damage and Tissue Repair
Numerous studies have provided evidence that elevated levels of fatty acids and of lipid species derived from fatty acids may cause cell damage resulting ultimately in cell loss by necrosis or apoptosis [41,42]. Enhanced damage and loss of cells elicits a regenerative response, which, in the liver, includes mitotic division of hepatocytes and in case of severe damage also the differentiation of hepatic stem cells. In the model, time-dependent changes of functionally intact cells in the i-th LU were described by the following equations: Here, η i is the fraction of functionally intact hepatocytes in the i-th LU, v i r and v i d denote the rates of cell regeneration and cell loss. v r was chosen to capture two opposing effects: The factor (1 − η) takes into account that the regenerative response initiated by immune cells rises with increasing cell damage, the factor η 2 takes into account that tissue regeneration requires interacting cells. The combination of the two effects entails that the rate of tissue regeneration v r is a non-monotone function with respect to η (represented by the blue curve in Figure 3A). The rate v d of cell damage follows a first-order kinetics with respect to η (represented by the black and red lines in Figure 3A), whereby the rate constant is the sum of a basal constant k d accounting for the continuous but small cell loss in the healthy liver, and an additional FFA-dependent term. The parameter γ determines the increase of v d if FFA exceeds the threshold FFA critical . The steady-state fraction of active hepatocytes is given by the condition v r = v d , which in Figure 3A corresponds to the intersection points between the green curve and the straight lines. If the damage rate exceeds a critical value, the fraction of active cells tends to zero (see Figure 3B). The numerical values of FFA critical and of the kinetic parameters k d , k r and γ (see legend of Figure 3) were chosen such that for the healthy liver, the life-span of hepatocytes amounted to 365 days (1 year), whereas for the cirrhotic liver and a fatty acid concentration that is twice as high as the critical concentration, i.e., FFA = 2·FFA critical , the life span is shortened to 26 days according to findings in Macdonald, 1961 [43].

Simulation of FFA-Induced NAFLD Progression
Simulations of FFA-induced NAFLD progression of the whole liver were performed by numerical solution of the coupled Equations (1)-(3) for each LU. As the progression of NAFLD takes place on the time scale of months and years, in these simulations the daily variations of lipid species (see Figure 2) were neglected and their concentration values were put to the quasi-stationary 24 h time averages. Functional heterogeneity of the liver was taken into account by randomly varying the kinetic parameters of the LUs by 20%, which corresponds to typical variations of blood flow [44] and tissue variations of protein abundances [45,46]. The simulations started at t = 0 with a healthy liver and a physiologically normal FFA release rate of v at = 1.17 mM/min to the plasma. At time t = 5 years, the release rate was increased to 1.76 mM/min. This resulted in a rise of plasma FFAs and the onset of liver steatosis, which was quantified as percentage of LUs having a TAG content larger than 30 mM. In concordance with lipidomics data of human NAFLD and NASH livers [47], the simulations yielded an about 6-fold initial rise of steatosis followed by a moderate decline with progressing NAFLD.
Severity of NAFLD at time t was evaluated by the total fraction of intact hepatocytes (TFH), At early time points after onset of the FFA challenge, LUs endowed with the lowest capacities for TAG storage and/or tissue repair are unable to establish a new balance between damage and repair and thus lose their functionality, i.e., η tends to zero. The metabolic failure of these units entails a reduced overall FFA-esterifying capacity of the liver, which causes an additional rise of the FFA plasma level and thus increases the FFA burden to the remaining intact LUs. This self-amplifying cycle of organ damage may stop if a remaining fraction of intact LUs exists that are endowed with sufficient metabolic and repair capacities to cope with the elevated levels of FFAs.
Although the magnitude of the random variations of model parameters was identical (standard deviation = 20%), the severity of NAFLD after 30 years displayed large differences (see Figure 4C). The reason for this heterogeneity is that the parameters determining blood flow rate, uptake and removal of FFA (k u , k β , k +TAG , k -TAG ) and damage and repair (k d , k r ) vary independently from one LU to another. Hence, the number of LUs endowed with the most favorable combination of these parameters (fast blood flow + high rate of FFA removal + low damage rate + high repair capacity) may differ from liver to liver. For the example in Figure 4A, the overall fraction of active hepatocytes dropped within five years to a moderately reduced and stable level of 80%. The clinical correlate of this type of NAFLD is a steatotic liver with low-grade inflammation. For the example in Figure 4B, the overall fraction of active hepatocytes dropped continuously to 45%. Clinically, this time course corresponds to disease progression through the stages steatosis → NASH → cirrhosis.
The example in Figure 4D illustrates how the liver may recover after cessation of the FFA challenge. Without lowering the release of FFAs into the plasma, the THF declines within 30 years to about 65%. This can be prevented by interrupting the FFA challenge, i.e., setting the release rate of FFA into the plasma back to the normal value (1.17 mM/min in the model). A short FFA challenge of several months, as it may occur during total parenteral nutrition [48] or transient insulin resistance due to infections [49], had virtually no impact on the hepatocyte fraction. The longer the FFA challenge persisted, the lower was the fraction of hepatocytes, which could be rescued from further damage. After five years of ongoing NAFLD, interruption of the FFA challenge restored only 77% of the functional intact liver mass. The time courses in Figure 4D are in good concordance with reported findings according to which the regeneration of fibrotic NAFLD livers after removal of the lipid load was generally incomplete and took many months [50,51].
In all simulations, steatosis was generally highest immediately after onset of the FFA challenge. The decrease of hepatic TAG with progressing NAFLD reflects the decrease of active cells capable of synthesizing and storing TAG. The more aggressive the NAFLD, the faster was the decline of steatosis.

Size and Length Scale of Intra-Hepatic Parameter Heterogeneity Influences NAFLD Progression
The simulations of NAFLD progression shown in Figure 4 were performed for livers composed of 100 LUs. With a liver volume of 1500 mL, this means an LU volume of 15 mL, which corresponds to a spherical LU with a diameter of about 3 cm. This is a typical length scale, at which significant spatial differences in liver steatosis are commonly assessed by means of MRI techniques [52,53]. In order to check how strong were the statistical differences between individual NAFLD time courses dependent on the size and length scale of parameter heterogeneity, we carried out simulations with livers composed of 10, 100 and 1000 LUs, corresponding to a characteristic length scale of 6.6 cm, 3 cm and 1.4 cm, respectively. As before, random parameter variances were put to 20%. As endpoint, we used TFH 30 , the total fraction of intact hepatocytes at 30 years after onset of the FFA challenge. For comparison, we carried out the simulations also for completely homogenous livers, where each liver was composed of identical LUs, the parameters of which were obtained by 20% random variation of the standard parameter set given in Table 1. Intriguingly, the homogeneous liver showed an "all-or-none" characteristics of NAFLD progression: 90% of livers had none or only marginal signs of NAFLD but the remaining 10% ended up in complete liver failure ( Figure 5A). In contrast, complete liver failure was not observed for the heterogeneous livers also composed of 100 LUs, but the distribution of TFH 30 values was much broader compared to the homogenous liver ( Figure 5C). Hence, intra-hepatic functional heterogeneity prevents complete liver failure but entails on average a more severe NAFLD progression as compared with a functionally homogeneous liver. Increasing the number of LUs, the distribution of THF 30 values became narrower and the mean THF 30 values were getting smaller than the THF 30 of the homogeneous liver (indicated by the red bars in Figure 5). From this, it can be concluded that the risk for NAFLD progression to severe forms increases with an increasing length scale of functional heterogeneity.

Patterning of Steatosis as Predictor of Prospective NAFLD Development
The heterogeneous distributions of intra-hepatic functional capacities result in different patterns of the regional TAG distribution during NAFLD progression. This raises the question of whether the pattern of regional steatosis in an earlier phase of NAFLD progression may allow a prognosis of prospective disease development. We addressed this question by monitoring the TAG distribution across the individual LUs at 5, 10 and 20 years after onset of NAFLD. Livers with a severe outcome (TFH 30 values lower than 55%) had on the average a strikingly higher proportion of non-steatotic LU s but the steatotic LUs had a higher mean TAG content than those in livers with a modest NAFLD outcome (TFH 30 values larger than 75%). These two features of TAG distribution discriminating between severe and modest outcome of NAFLD can be captured by a steatosis pattern score (SPS), SPS = (fraction LU s with TAG < 30) × (mean TAG in LU s with TAG > 30) (5) where, the first factor represents the fraction of LUs with a TAG content lower than 30 mM and the second factor is the average TAG content of the complementary fraction of LUs with TAG content larger than 30 mM. Plotting TFH 30 against the SPS ( Figure 6A-C) yielded significant negative correlations with increasing duration of the disease. For example, all livers with SPS < 16 at ten years after onset of NAFLD had THF 30 values larger than 60%, whereas all livers with SPS > 21 had THF 30 values smaller than 60%. Thus, the regional steatosis pattern at earlier time points of NAFLD progression appears to contain valuable information of the further progression of the disease. This conclusion holds without saying for conditions where the FFA challenge persists and other disease-promoting hits are absent.

Exacerbation of NAFLD by Additional Liver-Damaging Hits
In the preceding sections, simulated progression of NAFLD was exclusively driven by an FFA challenge without presence of further liver-damaging "hits." Such hits may be caused, for example, by changes of the gut microbiome, alcohol consumption, xenobiotica, drugs or viral infections. As proposed by the "second hit" hypothesis [54,55] and more recently by the "multiple hit" hypothesis [11], additional threats are necessary to drive an initially steatotic liver into NASH and even more severe stages of liver damage.
However, as suggested by the simulations shown in the preceding sections, for livers with an unfavorable intra-hepatic distribution of metabolic and repair capacities, such additional hits are not required to drive the liver into severe pathological states (see example in Figure 4B). Nevertheless, occurrence of additional hits may result in a dramatic acceleration of NAFLD progression. This is illustrated in Figure 7 for a liver that was exposed to three consecutive damaging hits with duration of about one year. In the absence of the FFA challenge, the hits resulted in a decrease of the total hepatocyte fraction to about 85%, which was reversible because the timely distance between the hits (~3 years in this example) was long enough to allow for a complete recovery. The FFA challenge of the same liver resulted in a partially damaged liver with a stable TFH 30 of about 63%. If the fatty liver was additionally exposed to the same hits as the healthy liver, TFH dropped gradually to eventually fall below 30%. Of note, the extent of the hit-induced drop of the hepatocyte fraction became smaller with increasing number of hits, a feature that was consistently observed in these simulations. Every hit removes a fraction of LUs with less favorable damage-counteracting capacities so that, on the average, the damage resistance and repair capacity of the surviving LUs was increasing.

Intra-Hepatic Distribution of Functional Capacities Influences NAFLD Progression
The most important finding of our computational study is that the heterogeneous intrahepatic distribution of metabolic and tissue-remodeling capacities is a strong determinant of NAFLD progression. Whether the onset of fatty-acid induced cellular damage in a fraction of liver regions may trigger a cascading damage depends on the functional capacities of the unaffected regions to effectively handle FFAs. We hypothesize that an unfavorable intra-hepatic distribution of functional capacities can be sufficient to drive simple steatosis to severe forms of NAFLD without occurrence of additional hits. However, as with all diseases, numerous additional risk factors can accelerate worsening of the disease. Our simulations suggest that even in livers with mild NAFLD, the impact of such additional "hits" is much more severe than in the healthy liver. Interestingly, the repeated occurrence of transient liver-damaging events represents a kind of selection pressure that favors the "survival" of LUs possessing the largest capacities to withdraw and repair tissue damage (see Figure 7D).

Cascading Failure as a Fundamental Principle of Disease Progression
Our study of NAFLD progression rests on the principle of cascading organ failure: loss of function due to the failure of parts of an organ/tissue has to be compensated by its functionally intact parts. This in turn exerts an additional load to the intact parts thus causing them to fail as well if adaptive mechanisms have reached their limit. In our model, the adaptive response was implicitly taken into account by the a priori statistical distribution of functional parameters across the LUs. It is reasonable to assume that the principle of cascading organ failure applies to the disease progression in other organs, e.g., heart failure due to long-lasting hemodynamic overload, liver or kidney failure caused by drug intoxication or insulin resistance and inflammation of adipose tissue caused by a TAG load exceeding the expansion limit.

Macro-Scale versus Micro-Scale Heterogeneity
Our mathematical model rests on the assumption that distinct spatial regions (that we named liver units-LUs) in the size range of centimeters exist, which differ from each other by blood perfusion rates, metabolic capacities of hepatocytes and the endowment with other cell types, in particular Kupffer cells and fibroblasts involved in defense and repair processes. This macro-scale heterogeneity has to be distinguished from the well-studied micro-scale heterogeneity among different zones of the acinus, commonly referred to as metabolic (functional) zonation. The existence of functional macro-scale heterogeneity of the liver can be inferred from the spatially variable distribution of blood flow rates, lipid accumulation and clearance rates of drugs and metabolites revealed by various imaging techniques. The molecular and cellular basis of these variations remains elusive. The vasculature of the liver may play an important role. For example, the pancreatic hormone insulin, which stimulates lipogenesis, has quite different concentrations in the various tributaries of the portal blood. Hence, hepatic territories low in insulin may explain areas spared by fat in a steatotic patient [56]. Concerning the contribution of regional differences in gene expression to intrahepatic functional heterogeneity, existing studies do not go beyond the single acinus although the technical prerequisites are now available [57]. Generally, it is well established that environmental fluctuations (extrinsic noise) affect the development of individual organisms and tissues [58]. Gene expression studies across mammalian organs suggest that random external events may largely influence the temporal trajectories of gene expression during organ development [59].

Spatial Patterning of Steatosis: An Additional Risk Factor?
Our model simulations suggest an association between the pattern of steatosis and the long-term outcome of NAFLD. A strongly contrasting steatosis patterning, i.e., the coexistence of non-steatotic and highly-steatotic regions, appears to be indicative for the later occurrence of severe NAFLD stages (see Figure 6). In the light of this finding, it could be worth to carry out a retrospective study relating clinical and histopathological data on NAFLD progression to the spatial pattern of steatosis revealed by ultrasound or nuclear magnetic resonance at the time point of diagnosis. If such study supports the modelderived hypothesis, spatially resolved patterning of steatosis by means of non-invasive imaging techniques might be included in the list of markers for the risk of a simple steatosis to progress to more severe stages of NAFLD.

Further Extensions of the Model
The presented mathematical model of NAFLD progression is a generic one in that fundamental processes involved in the initiation and progression of the disease have been included in the model as events that trigger each other. Gradually including the huge variety of intertwined molecular and cellular processes underlying these events may seem as a useful extension of the model. However, this extension should not be carried out according to the motto "put in everything that is known" but with a clear strategy that is oriented at a defined medical goal. One of the central medical goals in NAFLD research is to better understand and prevent the transition from simple steatosis to NASH. In this respect, the more detailed modeling of mechanisms and pathways involved in lipotoxic cell damage appears to be crucial in view of controversial ideas on useful pharmacological interventions. For example, abrogating the increase of hepatic lipogenesis (e.g., by inactivation of SREBP) has been advocated as a promising way to prevent steatosis and thus NASH [5]. According to our simulations (see Figure 2) and experimental findings [4], this might be the wrong way because enhanced lipogenesis is an important defense mechanism to lower the concentration of potentially toxic FFAs. On the other hand, excessive accumulation of TAG may induce endoplasmic reticulum stress and cell damage [60]. This Janus-faced effect of hepatic TAG synthesis on the progression of NAFLD calls for a more detailed and physiology-based model of the cellular lipid metabolism, which in particular includes the various processes involved in the synthesis, growth and degradation of lipid droplets [61]. Another important aspect of NAFLD progression, which was neglected in our model, is the upregulation of metabolic pathways to compensate for the decline of active hepatocyte mass. For example, it is well established that pro-inflammatory cytokines have a strong impact of on the lipid metabolism of the liver [62]. In the inflammatory phase of the disease, cytokine-stimulated upregulation of lipogenic pathways could facilitate the incorporation of potentially toxic free FAs into TAG.