Timescale analysis of a mathematical model of acetaminophen metabolism and toxicity Journal of Theoretical Biology

(cid:1) We have created a model, examining acetaminophen metabolism and related hepatotoxicity. (cid:1) We modeled multiple pathways associated with APAP metabolism. (cid:1) Using numerical, sensitivity and timescale analysis we have identi ﬁ ed key parameters. (cid:1) Analysis highlights a critical acet- aminophen dose in terms of the model parameters.. dose between safe and overdose cases, timescales for exhaustion and regeneration of important cofactors for acetaminophen metabolism and total toxin accumulation as a fraction of initial dose.

deaths associated with acetaminophen overdose in the U.S. almost doubled over a 4 year period, from 98 deaths in 1997 to 173 deaths in 2001 (Nourjah et al., 2006). In the UK, 90-155 people died per year between 2000 and 2008 with additional deaths due to acetaminophen being taken with other drugs (Hawton et al., 2011). This ease of availability and lack of awareness of its potential hazards means that acetaminophen is responsible for 80% of drugassociated cases of liver injury (Ostapowicz et al., 2002), and druginduced liver injury has become the most common cause of acute liver failure and subsequently transplantation in Western countries (Lee et al., 2003). Much of our understanding of the metabolism and toxicology of APAP comes from animal models, particularly rat and mouse. Interestingly there is considerable variation in toxicity between species (Davis et al., 1974).
APAP is taken orally and is absorbed into the blood stream. It arrives in the liver via the hepatic portal vein and moves through the liver mass to the central vein (Fig. 1). In this time, APAP is absorbed into the hepatocytes where it is metabolised. In the liver, hepatocyte function is determined by position relative to the portal vein, with functions differing if a hepatocyte is near the blood inlet (periportal) or outlet (centrilobular), an affect known as zonation and is present across all areas of the liver (Allen et al., 2005). APAP is metabolised in the liver primarily by the sulphation and glucuronidation pathways (Riches et al., 2009;Mutlib et al., 2006), while around 5% is metabolised, via oxidation, to form the toxic metabolite N-acetyl p-benzoquinone imine (NAPQI) (Duan et al., 2001). A detailed pathway diagram is shown in Fig. 2 and a simplified one used as the basis for the mathematical modelling is shown in Fig. 3. The sulphation pathway involves the conjugation of APAP with the cosubstrate 3'-phosphoadenosine 5'-phosphosulfate or PAPS. This cosubstrate is finite within the liver cell and at toxic doses we see PAPS levels fall (Sweeny and Reinke, 1988) and a saturation of the sulphation pathway, leading to higher metabolism through glucuronidation and oxidation. The cofactors associated with the glucuronidation pathway have a much higher capacity than those of the sulphation pathway (Reith et al., 2009) and we assumed in our modelling that the pathway does not saturate at clinically relevant, high APAP doses. Via the oxidation pathway, APAP is catalysed by select enzymes from a 'superfamily' of enzymes known as Cytochrome P450 (Patten et al., 1993). The main enzymes involved in this reaction in human cells are Cytochromes CYP2E1, CYP3A4 and CYP1A2 (Patten et al., 1993;Chen et al., 1998;Thummel et al., 1993), however, the sub-type and hence nomenclature of the enzymes varies by species when looking at animal models. Metabolism through oxidation produces NAPQI, a chemically reactive and toxic metabolite. NAPQI can be detoxified by GSH, an antioxidant which conjugates to NAPQI preventing binding with essential proteins and thus preventing damage to the liver. At sufficiently high doses, the sulphation cosubstrate, PAPS, can be exhausted, diverting quantitatively more APAP through the oxidation pathway, leading to higher amounts of NAPQI being produced. There are marked species differences in the sensitivity to APAP, e.g. rats are resistant to equivalent doses of APAP compared with humans, and this is due to a much greater capacity for sulphation and a lowered propensity for oxidation (Vaidyanathan and Walle, 2002). Oxidation has the effect of depleting GSH levels in the liver, through binding with NAPQI and hence greater levels of protein adducts are produced. GSH can also be depleted by individual factors such as alcoholism (Guerri and Grisolia, 1980) and anorexia (Kalsi etal., 2011) though this inter-patient variability is beyond the scope of the mathematical model to be presented in this paper.
It is broadly recognised that mathematical modelling now plays a significant part in the drug development process. A successful model provides a cost effective way of understanding and predicting drug efficacy and toxicology, thus offering a systematic means of guiding more focused, less exploratory, use of animal models. Despite acetaminophen being the subject of laboratory studies for many years, it is only recently that theoretical studies on the toxicology of paracetamol have been undertaken. One of Fig. 1. Structure of the liver (Frevert et al., 2005). Blood flows from the portal field (left) to the central vein. APAP in the blood diffuses into the hepatocytes and is metabolised.

Fig. 2.
A diagram of the cell scale metabolic network for APAP metabolism. The abbreviations are: APAP, acetaminophen; UGTs, UDP-glucuronosyltransferases; SULTs, sulfotransferase; NQO1, NADPH-quinoreductase; CYPs, cytochrome P450; APAP-G, acetaminophen glucuronide; APAP-S, acetaminophen sulphate; NAPQI, N-acetyl-p-benzoquinone imine; GSTs, glutathione S-transferase; GSH, glutathione; APAP-GSH, acetaminophen glutathione conjugate. Subscript 'B' denotes non-specific binding to a protein or lipid. Subscript 'P' denotes binding to non-specific protein (Diaz Ochoa et al., 2012). Blue boxes are non specifically bound products, yellow boxes are molecules, white boxes are isozymes, red boxes are protein bound molecules and green boxes are further metabolic systems not described in this diagram. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.) Fig. 3. Pathway Diagram for APAP Metabolism. APAP is metabolised through 3 main pathways, sulphation, glucuronidation and oxidation. CYP oxidation creates NAPQI, a harmful metabolite which can bind with essential cellular proteins within the hepatocytes if no GSH is present. Modelled species are APAP (P), NAPQI (N), PAPS (S), GSH (G) and Drug-Protein Adducts (C). the first mathematical models produced is by Reith et al. (2009), who focused on examining the kinetics of the glucuronidation and sulphation pathways using a 14 variable ordinary differential equation ODE) model and fitting to human data, specifically excreted products in the plasma. (Ochoa et al.) took a multiscale approach, combining a detailed cell based APAP metabolism model, comprised of 34 variables, with a whole body model to simulate actions in the liver and transport between organs. Both these models are rich in detail and parameter estimation, but their complexity prohibits investigation using more advanced mathematical techniques. Multi-compartmental models have also been tested by Ben-Shachar et al. (2012) who looked to create a model that would reproduce clinical and experimental data on APAP and metabolite levels in the plasma and urine. They looked to reproduce the data of (Prescott, 1980) examining APAP metabolism in human patients. Again, this model is complex and so it is difficult to apply mathematical analysis. Remien et al. (2012) investigate a simple model for APAP metabolism, utilising a tissue-scale model to predict biomarker levels, which can be used to estimate overdose amount, time elapsed since overdose, and likelihood of patient survival. In this paper we will present a cell-based model that describes the major pathways in the system, which is more detailed then the model proposed by Remien et al. but very much simpler than that of Reith et al. (2009) and Ochoa et al. This model will in fact be applicable to a broad range of drugs that are metabolised in the liver via (1) a non exhaustible pathway (i.e. glucuronidation), (2) an exhaustible pathway (i.e. sulphation) and (3) an oxidation pathway that leads to GSH binding and toxic conjugate formation. The resulting model is amenable to two forms of analysis. Firstly, to identify which parameters have the most affect on the predicted outcome through sensitivity analysis and, secondly, to derive relatively simple formula, using singular perturbation analysis, for factors such as critical initial dose and timescales for peak toxic activity. This will enable us to probe the model to gain great insight in to how individual mechanisms in the model can affect and influence these factors. Though the focus will be on APAP metabolism in humans, the modelling and analysis is applicable preclinical animal models also.
We seek to create a model that captures the most important aspects of APAP metabolism and toxicity at the cellular level. We then analyse the model both numerically and analytically in order to develop a better understanding of the interactions in the modelled system. We also wish to identify any data gaps which can then be pursued experimentally. In the next section we will derive the model. In Section 3 we present simulations, showing the metabolic responses to bolus doses of APAP and undertake parameter sensitivity analysis. In Section 4 we perform a detailed timescale analysis, to derive formula characterising APAP metabolism. Finally we summarise the key results and discuss future work in Section 5.

Model background
We focus on the metabolism of paracetamol within a single hepatocyte, aiming to capture the main dynamics of APAP metabolism while maintaining enough simplicity that analytical progress is possible. The full metabolic process is summarised in Fig. 2 and, as stated before, broadly separates into three pathways. Describing all of the pathways illustrated in Fig. 2 would lead to an extremely complex model involving 20þ state variables and many more parameters. Instead, as a first approximation, we bundle all the pathways in the glucuronidation route into a single pathway and likewise for sulphation and oxidation. The reduced pathway diagram used for the model is shown in Fig. 3. We assume for sulphation and glucuronidation that the first reaction down each pathway is non-, or negligibly, reversible, so that events downstream do not directly affect paracetamol metabolism. For the oxidation pathway, we assume a single generic CYP is involved which represents the combined actions of CYP2E1, CYP3A4 and CYP1A2.

Model description
We use mass action laws to derive a system of ordinary differential equations that describe the dynamics over time of the different pathways illustrated in Fig. 3. The resulting model is the same as that presented, but not studied, in Williams et al. (2013), we will nevertheless outline the model derivation. The model variables are listed in Table 1 and we note they represent quantities per cell.
Our model assumes an initial bolus dose being delivered. The metabolism depends on the size of the initial dose. At regular doses the majority of APAP will be metabolised by the sulphation and glucuronidation pathways (Duan et al., 2001). APAP (P) undergoes sulphation by reacting with the PAPS enzyme (S) at rate k S SP, where k S is the rate constant associated with the metabolism of APAP by PAPS, ultimately forming APAP-S. In humans, PAPS is exhaustible and so at high doses may, in some situations, see a saturation of the pathway. We define the rate constant for the production of PAPS by the liver as b S and the rate constant for the natural decay as d S . In contrast, we assume that the enzymes involved in glucuronidation are not exhaustible and are present at an approximately constant concentration, hence APAP metabolism along this pathway is in affect a natural decay at rate k GP .
The remaining APAP is metabolised via the oxidation pathway creating NAPQI (N). We assume that cytochrome P450 enzymes are present continuously at an approximately fixed concentration, so that the oxidative pathway is described as a further "natural decay" term, k 450 P. This reaction is assumed reversible at rate k NN .
NAPQI is assumed to be metabolised via one of two pathways. The first is by reaction with the antioxidant GSH (G) at a rate k GSH NG. At normal doses of APAP we expect to see nearly all of the NAPQI produced being detoxified by this pathway. Conjugation with GSH renders NAPQI harmless and it is excreted from the body with no ill effects. In our model GSH is assumed to be constitutively produced at a constant rate b G and naturally decays at rate d GG . In fact, the production and regulation of GSH production is quite complex, being released from skeletal muscle (Bilinsky et al., 2015) and regulated as an adaptive mechanism by NRF2 (Klaassen and Reisman, 2010); at the level of detail of the current model we assume that constant b G is a reasonable starting point for modelling single doses. The second pathway has NAPQI creating drug-protein adducts (C) at a rate k PSH N. This binding to cellular macromolecules can result in cell death if the proteins that are bound are essential for cell function/viability. We do not consider the downstream events caused by drug-protein adducts and the variable C represents the total accumulated amount of a toxic reaction (we therefore hereon refer to C as toxins in that they are capable of inducing cell death). We arrive at the following model describing the pathways in Fig. 3 and including the stated assumptions: We assume in this study that the drug is introduced into cells as a single bolus dose at t ¼ 0 at a concentration P S . The cells at this point are assumed to be at pretreatment steady-state level. The initial conditions for this system are thus Table 2 lists the model parameters and their estimated values for the standard simulation. Where possible, we obtained their values from the literature and any remaining parameters through repeated simulation, so that the numerical results matched reasonably well with similar simulations from Remien et al. (2012). It is generally considered that anything more than 4 g taken at once is considered an overdose, so we use 4 g as our safe dose case (Craig et al., 2012) (though it is recommended to take no more than a 1 g dose at 4 h intervals).

Simulation
Our aim is to understand the effect of dose on both NAPQI production and timescales of events in APAP metabolism. We solve the system of Eqs. (1)-(5) using the MATLAB routine ode15s, a variable order backwards difference method. Unless otherwise stated, we use the parameter values listed in Table 2.
We first examine the single 4 g dose case, i.e. a daily dose in a single bolus. We expect GSH levels to remain non-negligible to ensure a safe low-level conjugation of NAPQI. Consequently, protein adducts will then stay at very low levels. Both of these features can be observed from the simulation in Fig. 4 (left column).
It can be seen that neither GSH or Sulphation levels drop to zero, indicating that all APAP in the system is being dealt with effectively. We do see a rise in NAPQI, however overall levels are extremely low relative to our overdose case and therefore do not pose any great risk. The same can be observed for the protein adducts, which remain at low levels compared to the overdose case.
For the overdose case of 16 g, a likely outcome is that both GSH and sulphate levels will become exhausted at some stage of the metabolism process. This indeed occurs as can be seen in Fig. 4 (right column). Sulphates drop very rapidly to a near zero level and take a long time to recover; this means that proportionally more APAP will be conjugated into NAPQI. This leads to a rapid drop in GSH to negligible levels that are sustained for a period of about 40 h. This rise in NAPQI and subsequent depletion of GSH results in a high level of formation of protein adducts in comparison with our safe dose simulation. We note that a 4 Â increase in dose leads to an almost 10 4 Â increase in accumulated protein adducts. Fig. 5 shows the affect of the initial dose on the total amounts of toxic protein adducts produced, presented as C 1 =P S , where C 1 represents the steady state level i.e. C-C 1 as t-1 and the ratio C 1 =P S represents the fraction of adduct molecules produced per APAP molecule. At levels just slightly above a safe dose of 4 g it can be seen that the amount of protein adducts in the system rises rapidly. This rapid increase in protein conjugate formation displays how dangerous overdoses involving APAP are. Small increases in the dose above what is considered "safe" lead to huge increases in the protein adducts being produced, which in turn can lead to extensive damage to the liver. This threshold behaviour is due entirely to the level of GSH depletion of which leads to the fraction of protein adducts produced increasing 1000-fold over a 3-5 g dose (we note that in Remien et al. the lowest doses for patients receiving treatment is about 5 g). The sensitivity of the model solutions to parameter change is explored in the next section, whilst the key parameters governing the threshold dose are established in the analysis of Section 4. We note that as P S -1, the sulphation pathway becomes less significant and it follows that Simulations investigating the effect of smaller regular doses are shown in Fig. 6, in particular those in the left column represent a typically prescribed 1 g dose at 4 separate 4 h intervals over a 5 day period. Here, we observe NAPQI progressively building up in the initial days before settling to a periodic profile. Protein adducts increase linearly, although total levels still remain negligible.
The right hand side of Fig. 6 plots a higher than recommended chronic dose case, this time with the patient taking 1.5 g of APAP every 4 h. This increase in APAP leads to a rapid depletion of GSH resulting in NAPQI and conjugate levels two orders of magnitude higher than in the 1 g case. NAPQI and protein adducts both rise rapidly (after a day) due to the lack of GSH in the system to safely deal with the NAPQI present. The plots once again show dramatic increase in toxic effects (represented by an increase in adducts, C) following a modest overdose.

Parameter sensitivity analysis
The results in Section 3 demonstrated a notable sensitivity to dose. In this section we seek to establish the sensitivity of the See (5) (1) 4 g dosage, standard single dose assuming 80% of dose reaches liver.
(2) Assuming initial PAPS is 10% of standard APAP dose i.e. bS dS ¼ P0 10 , and initially sulphation and glucuronidation are about the same, i.e. k S ¼ kG dS bS i.e. amounts to 47.5% of APAP processing initially.
(3) Equal to kG 9:5 i.e. we assumed only 5% of APAP is oxidised initially. (4) Assumed k N ¼ k450 10 i.e. forward reaction is dominant. (5) Assuming at normal GSH concentration, bG dG , only 1% of NAPQI binds with the hepatocytes, i.e. k PSH ¼ 0:01 kGSH bG dG . Parameters marked with [*] indicate parameters chosen by us to produce physiologically realistic results. model solution to changes in parameter values. To do this systematically we used the Latin Hypercube method implemented using the "lhsdesign" routine in MATLAB. To produce the results that followed, the routine was set up to run 500 iterations, which randomly selects parameters between set limits of 3 Â and 1 3 Â their original value. We used, for the sensitivity test, the total  Fig. 4. Plot of the evolution of, from top to bottom, APAP, PAPS, NAPQI, GSH and protein adducts respectively. The units in each graph are mol/cell, noting the two orders of magnitude difference between the levels in N and C. Here 4 g (left) and 16 g (right) correspond to P 0 ¼ 1:32 Â 10 À 13 and P 0 ¼ 5:28 Â 10 À 13 mol/cell, respectively. accumulated protein adducts C 1 (i.e. CðtÞ as t-1), where we plotted this against each of the model parameters. We look for trends in the resulting graphs, indicating higher or lower numbers of protein adducts in response to a change in parameters. To confirm our observation we also examined the Sobol indices to estimate the sensitivity of variance of the model output, C, to the variance of the parameters (Saltelli, 2008). Defining indices S i (the first order effect) and S T i (the total effect index) to be the conditional expectation divided by the unconditional variance and the total output variation due to a given parameter respectively. Then S T i À S i ¼ 0 indicates that a parameter has no affect on the variance of the model output.
Shown in Fig. 7 are the results of the sensitivity analysis for the safe dose of 4 g. We observe that most of the graphs do not show any sort of trend in response to differing parameter values except that the k 450 (oxidation rate constant) graph shows an obvious upward trend in protein adducts whilst a downward trend is observed for k G (glucuronidation rate constant). Here, the Sobol indices are found to be S T k G À S k G C 0:35; S T k 450 À S k 450 C0:32 and S T b G ÀS b G C 0:2, while all other values are less than 0.05 confirming the visual analysis of the parameter sensitivity. Interestingly, the indicated sensitivity to b G is not present in the 16 g case, suggesting that this is likely to be an important parameter when doses are near to the "critical level". However, as the Sobol Indices for b G in the overdose case indicate that it has no significant affect on the model output, we do not feel that any further analysis is necessary for this parameter. The sensitivity analysis for the overdose case of 16 g is shown in Fig. 8. Again we see that changes in the value of k 450 and k G produce the most distinct trends in the model.
In the overdose case, S T k G À S k G C 0:16 and S T k 450 À S k 450 C 0:06, while all other values are again less than 0.05. These are the only 2 parameters with a notable affect on the model outcome in the overdose regime.
This analysis suggests that the key mechanisms that govern paracetamol metabolism are glucuronidation and oxidation, where increasing k G reduces toxicity and increasing k 450 enhances it. In the parameter range investigated, PAPS contributes only up to about 10% of APAP metabolism, whereby sulphation is a secondary process in humans; we note that the sulphation for rats lies outside the parameter range investigated. Figs. 9 and 10 show the dependence of the total toxins produced, C 1 , on the two most sensitive parameters k 450 and k G , for the safe and overdose cases. The results were generated from running the simulation to approximate C 1 ðt ¼ 50Þ, we found this length of time sufficient to reach a steady state. From Fig. 9 we see that increased k 450 will lead to more APAP being oxidised instead of being metabolised by sulphation or glucuronidation. This will cause a rise in the amount of NAPQI in the system, putting more strain on the GSH pathway. We anticipate that a higher value for k 450 will lead to more protein adducts being present in the system and therefore increase the risk of liver damage.
The safe dose response shows a steady increase in conjugate levels initially, followed by a rapid rise in conjugate levels being produced with total protein adduct formation increasing by over one order of magnitude. A less dramatic rise in protein conjugates is observed for higher k 450 values. For our overdose case we see a much faster rise in the total protein adduct formation in response to higher k 450 levels. We see an increase of approximately three orders of magnitude in response to higher values of k 450 . After the initial rapid increase in total protein adduct formation, higher values of k 450 have a much lower affect on C 1 . Once GSH is depleted in our system, all NAPQI that is oxidised will produce protein adducts, the rate at which these protein adducts can be formed is then dependent on how quickly NAPQI can be oxidised, this rate is k 450 . This suggests that after GSH is depleted fully, conjugate production will be proportional to k 450 . The rate of APAP to NAPQI metabolism can be affected by other factors such as caffeine consumption (Carter, 1987;Lee et al., 1996) and, for example, consumption of anti-convulsant drugs (Bray et al., 1992) which would result in a higher value of k 450 .
In Fig. 10 we observe, as expected, a decline in toxins produced as k G increases. As with k 450 , there is a fairly sharp transition between high and low toxicity at a certain value of k G . We note that a 10-fold increase in k G is required in the overdose case ðk G $ 18=dayÞ to produce minimal toxic levels in comparison to the safe doses ðk G $ 1:5=dayÞ. Furthermore, for k G $ 0/day there is a 10-fold difference in C 1 levels. This is due simply to more APAP being present for longer in the overdose case. The critical role of GSH exhaustion is highlighted in Fig. 11, which plots the numerically predicted minimum value against parameters k G and k 450 . Of particular note is how the value of k G and k 450 at which the sharp jumps occur correspond to jumps in Figs. 9 and 10. Fig. 11 plots the minimum GSH levels in the cell against k G and k 450 . As k G increases we see a rise in the minimum GSH level of 2 orders of magnitude. This suggests that if the glucuorindation rate drops then GSH could fall low enough to allow protein adducts to form, if for example a person has a genetic or environmental deficiency e.g. co-medication, that reduces the amount of glucuornidation cofactor it could be dangerous for them to take paracetamol, even in safe doses. We observe in the overdose case that only very large values of k G have a non-negligible affect on minimum GSH levels. An increase in k 450 leads to a drop in GSH of over 3 orders of magnitude in the safe dose case, at the normal value of k 450 ¼ 0:315=day, minimum GSH levels remain high in the cell. However, an increase in k 450 leads to lower GSH levels which could lead to the formation of protein adducts. Therefore, increased k 450 can potentially lead to liver damage via protein adduct formation even in safe dose cases.

Cellular dose variation
The structure of the liver lobule means that cells closer to the portal vein are likely to receive more of the drug. As a consequence, there will be a distribution of drug dosage between cells in the liver. Some cells which receive higher doses are more likely to be damaged than others. Furthermore, differences in microenvironment due to proximity to blood vessels and oxygen gradients could also affect drug metabolism. The effects of the micro-environment will be subject to a future publication, and here we investigate how the spread of drug dosing effects the probability of cell death given a dose.    Table 2. In Fig. 12, we assume that a dose of paracetamol is log-normally distributed between the hepatocytes in the liver. We assumed that a lethal dose for cells (p L )is 5 times the daily safe dose (Remien et al., 2012) and we plotted the probability p 4 p L , given a mean dose log ðp s Þ and variance σ 2 , against mean dose. We observe that higher standard deviations lead to a less sharp profile. It is expected that 70% total cell death will lead to the death of the patient (Remien et al., 2012), in our simulations we see that this occurs from $ 7 Â 10 À 13 mol=cell (approx. 5 times the standard dose) to $ 9 Â 10 À 13 mol/cell (approx. 7 times the standard dose). Interestingly, we also see that greater variation between hepatocytes leads to more deaths at lower doses, but less death at higher doses. This suggests that variation is a positive property for the population on average, for survival against a very large, single dose. However, this does not necessarily mean it is a positive property for the individual.

Timescale analysis
In the previous section we were able to get some insight into how certain parameters effect the predicted toxicological outcome. In this section we will employ singular perturbation theory to get a much better analytical understanding of APAP metabolism according to the model. Close examination of Fig. 4 reveals the existence of distinct timescales, starting with a rapid decline in PAPS and GSH followed by longer timescales for recovery. To apply this theory we first non-dimensionalise the system of Eqs. (1)-(5). Using the data values in Table 2, we express the new parameters in terms of a single small parameter ϵ (i.e. ϵ⪡1), which will be exploited in the analysis. We will summarise the main results here, and we refer the reader to the supplemental material for full details.

Non-dimensionalisation
To aid the analysis we rescale our variables in order to eliminate units, which allows comparison of variables and parameters in terms of their magnitude, so that the dominant and negligible mechanisms can be systematically identified. Since glucuronidation is the dominant metabolism route for APAP, we rescale time Fig. 9. Total protein adduct formation against k 450 for the safe (4 g, dashed) and overdose (16 g, solid) cases. The dotted line indicates the standard value corresponding to data in Table 2. Fig. 10. Total protein adduct formation vs. k G for the safe (4 g, dashed) and overdose (16 g, solid) cases. The dotted line indicates the standard value of k G found in Table 2.  Table 2. with parameter k G ; using the value in Table 2, the dimensionless timet ¼ 1 thus represents about 8 h. We rescale PAPS and GSH with their untreated levels and rescale APAP, NAPQI and protein adducts to a reference value P 0 which represents the liver cell level of a 4 g dose i.e. P 0 ¼ 1:32 Â 10 À 13 mol/cell. The rescalings are thus and we note the standard dose concentration P 0 corresponds tô p ¼ 1. The dimensionless system of equations is then dp dt ¼ Àα where the rescaled parameters are listed in Table 3. The third column of Table 3 lists the value of the parameter, and for the purpose of the analysis we will rewrite them in terms of the small parameter ϵ ¼ k 450 =k G C0:1 guided by magnitudes indicated in the 4th column; thus starred values in Eqs. (7)-(11) are defined aŝ k N ¼ ϵ 2k n N andα S ¼α n S . These dimensionless variables are subject to the initial conditionŝ pð0Þ ¼ P S ;ŝð0Þ ¼ 1;nð0Þ ¼ 0;ĝð0Þ ¼ 1;ĉð0Þ ¼ 0; recalling that P S ¼ 1 represents the 4 g dose case. Henceforth, we will drop the hats and the n's for clarity. In Section 4.2 we provide an overview of the main mathematical results and we then give biological interpretations in Section 4.3.

Application of singular perturbation theory
The system (7)-(11) will be analysed in the limit ϵ-0. Using singular perturbation theory we can perform this analysis systematically and formally reduce the full system to a sequence of timescales in which the system reduces to a simpler solvable one in each timescale. This will enable us to identify when a particular process is important and determine an approximation to key quantities such as critical dose in terms of the model parameters. Full details of the analysis is given in the supplementary material and we present only the "highlights" below. A summary of this analysis and the important timescales and events is provided in Section 4.3. We note that toxic levels of protein adducts will be c ¼ Oðϵ 2 Þ as shown in Fig. 5.

t ¼ Oðϵ 3 Þ
On introduction of the APAP bolus there is a rapid adjustment over t ¼ Oðϵ 3 Þ, the first 30 seconds or so, in which NAPQI is produced at very low levels. Denoting variables in this timescale with a superscript n, we write t ¼ ϵ 3 τ n ; p ¼ p n ; s ¼ s n ; n ¼ ϵ 4 n n ; g ¼ g n ; c ¼ ϵ 6 c n : These rescalings are then substituted into our dimensionless equations (7)-(11), subject to p n ¼ P S ; s n ¼ 1; g n ¼ 1; n n ¼ 0 and c n ¼ 0 at t n ¼ 0. In each timescale we seek solutions of the form: pðτ n Þ ¼ p n 0 ðτ n Þþϵp n 1 ðτ n Þþϵ 2 p n 2 ðτ n Þþ⋯ and likewise for the other variables. Substituting these expansions into our equations we obtain to leading order p n $ P S ; s n $ 1 and g n $ 1 (correction terms can be found in our supplementary material) and In this short initial timescale, APAP, PAPS and GSH remain relatively unchanged and NAPQI equilibrates to a negligible Oðϵ 4 Þ level. As t n -1, NAPQI settles to n $ ϵ 4 P S =α G À Á and c $ ϵ 6 k PSH P S τ n =α G À Á . We note here that as τ n -1; n $ ϵ 4 P S =α G , this represents the amount of NAPQI formed if PAPS and GSH remain at their pretreatment levels. There is no change at leading order of p; s and g, however the correction terms become Oð1Þ at τ n ¼ Oðϵ À 2 Þ i.e. at t ¼ OðϵÞ.

t ¼ OðϵÞ
It is on this timescale at which sulphation is most prominent.
We introduce t ¼ ϵτ and the relevant rescalings are  Substituting the expansions above, p $ p 0 þ ϵp 1 etc. into (7)-(11) and solving leads to p $ P S þϵ 1 ϕ S ðe À α S ϕ S P S τ À 1ÞÀP S τ ; s $ e À α S ϕ S P S τ ; g $ 1 þ ϵðÀϕ G P S τÞ; In this timescale, we see that sulphate levels drop rapidly whilst APAP is relatively steady. Biologically this is due to the conjugation of APAP and PAPS, leading to declining PAPS levels in the cell. The parameters used suggest that the pretreated PAPS concentration is OðϵP S Þ so, at best, sulphates are only able to metabolise an OðϵÞ fraction of the drug. There is also an increase in protein adducts, although they are still only present in very low amounts. We note as τ-1; p $ P S À ϵ ϕ À 1 S þ P S τ , where ϵ=ϕ S represents the amount of APAP being metabolised by the sulphation pathway. There is a transition timescale t ¼ ϵ η 1 ðϵÞþOðϵÞ, where η 1 ¼ lnð1=ϵÞ=α S ϕ S P S , in which sulphate reaches a minimum constant level, namely s $ ϵδ S =α S ϕ S P S ; sulphation makes no further contribution to APAP metabolism at leading order. The expansion breaks down when τ ¼ Oð1=ϵÞ, corresponding to t ¼ Oð1Þ, when APAP concentration starts to significantly drop.

t ¼ Oð1Þ
In this timescale, we have two separate divergent cases. One in which we have sufficient amounts of GSH in the system to conjugate NAPQI, the other is characterised by a rapid drop in GSH and potential toxin build up. The critical dose at which the two cases diverge is such that, P S oP n S can be classified as "safe" and P S 4 P n S can be considered a potential overdose. We note here that we have assumed that δ G a 1, we will omit details for the coincidental case of δ G ¼ 1 (i.e. δ G ¼ k G in dimensional terms).
In both cases, we adopt the following rescaling: t ¼τ; p ¼p; s ¼ ϵs; n ¼ ϵ 4ñ ; g ¼g ; c ¼ ϵ 3c : Expanding these variables in the usual way, and solving the resulting system yields p $ P S e Àτ Àϵ e Àτ ð c $ k PSH Zτ 0ñ ðτÞ dτ: Here, APAP is metabolised such that p $ P S e À τ (due to glucuronidation at leading order) and that PAPS is recovering, noting that s ¼ OðϵÞ and therefore is not contributing to APAP metabolism at leading order. We also note thatc is unsolvable in this timescale, but we can deduce behaviour asτ-τ n (see below), as explained in the supplementary material.
The divergence depends on the function Ψ ðτÞ ¼ 1 þ ϕ G P S δ G À 1 ðe À δ Gτ À e Àτ Þ; 8τ 4 0 whereby if Ψ ðτÞ 4 0, 8τ 40, thenñ andg remain positive and Oð1Þ, this is our safe dose case. If atτ ¼τ n , such that Ψ ðτ n Þ ¼ 0, thenñ-1 in finite timeτ-τ n whilstg-0. The divergence condition ðP S ¼ P n S Þ is determined by assuming that Ψ ðτ n Þ ¼ 0 is a turning point atτ ¼τ n , i.e. solving Ψ ðτ n Þ ¼ 0 and Ψ 0 ðτ n Þ ¼ 0 simultaneously leading to t n ¼ lnδ=ðδ À 1Þ. We note that the safe and overdose cases can be connected smoothly by analysis in the region of P S ¼ P n S þ ϵθ, where θ C Oð1Þ. The results are omitted as they are not of biological significance other than it reveals that the jump region observed in Fig. 5 is of OðϵÞ ¼ Oðk 450 =k G Þ in size.
In the overdose case, when P S 4 P n S , breakdown occurs when t $ μ 1 ðϵÞ, where μ 1 ðϵÞ is defined such that Ψ ðμ 1 ð0ÞÞ ¼ 0; andg ¼ OðϵÞ andñ ¼ Oð1=ϵÞ. Here, μ 1 ðϵÞ is the time at which hepatocytes no longer have an effective means of dealing with NAPQI. It is straightforward to show that μ 1 ðϵÞ is a decreasing function of P S and d g , i.e. more drug and less glutathione reduces the time interval, as expected. We further note that given Ψ ðμ 1 Þ ¼ 0 and Ψ 0 ðμ 1 Þ o 0 we can show that ϕ G P S e À μ 1 4 δ G ; this result is utilised in Section 4.2.5. In the overdose case, breakdown occurs when t $τ n ¼ μ 1 ðϵÞ, where Ψ μ 1 ð0Þ À Á ¼ 0 and Ψ ðτÞ ¼ OðϵÞ, so thatg ¼ Oð ϵÞ andñ ¼ Oð1=ϵÞ, this is discussed in Section 4.2.5.

Safe dose case ðP S o P n
S Þ Here, the drug decays exponentially (predominantly by glucuronidation) andg ¼ Oð1Þ throughout, i.e. GSH is able to handle the NAPQI being produced. Meanwhile sulphate cofactors are recovering but only at very low levels. Protein adducts attain their maximum level i.e. Oðϵ 3 Þ, namely There is a further timescale at t ¼ lnð1=ϵÞþOð1Þ in which the sulphation factors, now Oð1Þ, continue to recover, and return to pre-treatment state.

t ¼ μ 1 ðϵÞþOð1Þ (overdose)
GSH and NAPQI continue to drop and rise, respectively, over a series of intermediate timescales until the current one describing the time period at which GSH is at its minimum level. We rescale our variables as follows: We then expand our variables as before, substitute them into (7)-(11) and solve to find p $ P S e ð À μ 1 À τÞ ; In this timescale APAP levels continue to drop while sulphates remain steady. Protein adducts approach their maximum level while NAPQI production begins to slow and GSH levels begin to rise as APAP levels decline. The solutions in this timescale breakdown as τ-μ 2 ðϵÞ À , with μ 2 ð0Þ ¼ lnðϕ G P S =δ G ÞÀμ 1 ð0Þ, where and n ¼ OðϵÞ. After this timescale, NAPQI levels begin to decline. As τ-μ À 2 , c attains its maximum value to leading order, i.e. c 1 $ ϵðP S e À μ 1 ð1 Àe À μ 2 ÞÀμ 2 δ G =ϕ G Þ. We can show that the amount of protein adducts increases with P S (i.e. higher initial dose) and ϕ G (less GSH present) as would be expected.

t ¼ μ 1 ðϵÞþμ 2 ðϵÞþOð1Þ(overdose)
This timescale follows a series of intermediate timescales in which GSH rapidly recovers and NAPQI diminishes. Here, the rescalings are and proceeding as before Here we see that APAP levels continue to drop exponentially, allowing PAPS levels to rise exponentially. GSH levels are now Oð1Þ and will soon recover to its pretreated level, whilst the tiny amounts of NAPQI that remain rapidly decrease. We now have GSH returning to pretreated levels as NAPQI diminishes.
After this, the only timescale of significance is τ o ¼ lnð1=ϵÞþOð1Þ, whereby p-OðϵÞ and s-Oð1Þ, i.e. their pretreated levels. Fig. 13 shows the evolution of the dimensionless APAP, PAPS and GSH concentrations against dimensionless time in an overdose case (left). As expected, the agreement improves as ϵ decreases (right).

Timescale analysis summary
Here we summarise the important events and timescales from the previous section, expressing the key dimensionless quantities in their dimensional form.

Critical paracetamol concentration
In Section 4.2.3, where t ¼ Oð1Þ we observe a divergence between our safe and overdose cases. This divergence occurs at a critical concentration where P n S ¼ 1:47 Â 10 À 13 mol=cell using the data available in Table 2. We note 4 g translates to a concentration of 1:32 Â 10 À 13 mol=cell and our divergence happens at a point 11% above this dose. This highlights the relatively low tolerance the liver has in response to large bolus doses of paracetamol.

Exhaustion of sulphate
Our analysis shows that sulphate is exhausted in the intermediate timescale between 4.2.2 and 4.2.3. The approximate timescale for exhaustion of sulphate is t $ k 450 k G k S P S lnðk G =k 450 Þ; which using the data is t $ 12 min for a 4 g dose. After this point the pathway saturates and we a greater proportion of APAP being metabolised into NAPQI, impacting GSH levels. We note that the estimate is only logarithmically accurate and will not be as precise as those in Sections 4.3.1, 4.3.4, and 4.3.5 are; nevertheless it makes explicit how much faster PAPS is exhausted in response to an increased drug dose.

Sulphate recovery
In both safe and overdose cases, we see sulphate recover at t ¼ lnð1=ϵÞ, in dimensional parameters this is Using the data this equates to about 40 h after ingestion for a 4 g dose; though we note, like that of Section 4.3.2, this estimate is only logarithmically accurate. Sulphate recovery is a long term process and the liver takes a long time to recover from a high paracetamol dose. In the case where a person uses paracetamol chronically to deal with pain then this long recovery time could impact how well the liver can deal with multiple doses. We note, as expected, that the recovery time is extended with dose, but in a sublinear fashion.

GSH depletion
In our overdose case, when P S 4 P n S we observe a collapse in GSH levels at t $ μ 1 (Section 4.2.3). Where μ 1 satisfies the implicit This equation cannot be solved directly to find μ 1 but given values of the parameters, the equation can be solved using the Newton-Raphson method. Using the data in Table 2 gives μ 1 % 0:046 for the overdose case, which using dimensionless parameters is μ 1 % 0:138 which is in good agreement with the numeric values shown in Fig. 13. We then show that mathematically we can improve our estimate by reducing the size of ϵ. Similarly, we find t $ μ 2 % 2:358 in the overdose case, again providing good agreement with the numeric values shown in Fig. 13. In terms of dimensional parameters this gives us μ 2 % 0:79, which is discussed further in Section 4.3.5.

GSH recovery
Again looking at the case where P S 4 P n S , the time for GSH recovery is given by which is approximately 8.9 h for a 4 g dose and 20 h for a 16 g dose. People regularly taking high doses of APAP can cause damage by not allowing time for GSH recovery and subsequently protein adduct formation could be high. Again, the plot in Fig. 13 shows how this estimate of GSH recovery is accurate for our model and how smaller values of ϵ (i.e. a decreasing k 450 =k G ratio) increase the accuracy of our estimate.
4.3.6. Total protein conjugate formation, C 1 The total concentration of drug-protein adducts in the overdose case, P S 4 P n S , is We note that the accumulated drug-protein adducts total is unaffected (to leading order) by parameters associated with PAPS. Moreover, we can show that C 1 increases with an increasing initial APAP dose and CYP reaction rate, and decreases in response to an increasing GSH production rate and glucuronidation reaction, as expected. Fig. 13 shows that this offers a fair prediction of maximum drug-protein adduct levels. We note that no parameters associated with sulphation have an affect on the final accumulation of protein conjugate formation and suggest that the sulphation pathway is unlikely to be a suitable target for an effective new treatment against the toxicological effects of an APAP overdose. However, we should note that even though sulphates are "exhausted" by time t $ lnðk G =k 450 Þ=k G P S it is still removing APAP at around the same rate as the oxidative pathway between timescales 4-10 (see supplemental material).

Discussion
In this paper, we have derived a cell scale mathematical model which describes the metabolism of APAP in hepatocytes. In order to obtain insights into this system using analytical methods, we simplified the full metabolic pathway to one that still maintains the three major pathways.
The simulations demonstrated that the model captures the expected dynamics of metabolism and, in particular, the distinguishing dynamics between the safe and overdose cases. We observe the expected drops in both sulphate and GSH levels in the safe dose case and our overdose simulations have both pathways dropping rapidly to very low levels, which is what we expect from clinical observation. The results show that a 4 Â dose of APAP can lead to a 100-10,000 Â increase in the amount of protein adducts being formed.
Our sensitivity analysis has enabled us to identify the most sensitive parameters in our model, we can use these to guide the research of biologists which will then provide further insight and help us to refine our model. The analysis in Section 3.2 showed that the key parameters are k G (the rate constant for glucuronidation) and k 450 (the rate of oxidation); the other parameters have secondary effects on the dynamics and, in particular, the sulphation pathway is less influential than glucuronidation and oxidation. There is ongoing work by the authors examining adaptive responses to chronic dosing. For example, if certain pathways become up-regulated in response to mild liver stress caused by APAP, then these sensitive parameters may be one of the contributing factors.
It can be seen that system operates over a number of distinct timescales. At t ¼ OðϵÞ $ 45 min we see sulphate levels begin to decline in response to the APAP present. As time progresses, we observe that sulphates begin to decline, and by t ¼ Oð1Þ $ 8 h we see that sulphates have become exhausted as they drop by an order of magnitude (i.e. S ¼ OðϵÞ). At this stage in our analysis, we see a critical divergence between the safe dose and overdose cases at an initial paracetamol dose of 4.54 g (using data from Table 2). We also are able to identify timescales for exhaustion and recovery of GSH and sulphates (details of which are available in the supplementary material). Of course there can be considerable individual variability that can affect the critical dose level. The sensitivity analysis has enabled us to deduce that changes in k G and k 450 have the largest impact on the dynamics of the system. Further to this, the asymptotic analysis of Section 4 has allowed us to express key quantities (critical concentrations and timescales) in terms of relatively simple formula (Section 4.3), so the effect of varying parameters can be explicitly observed. Such methods have broad application and are somewhat underused in the study of mathematical models in pharmacology.
Our parameter selection is good but there are gaps in the current literature that highlight a need for more data on the metabolism of APAP in humans. While literature is available which has allowed us to begin parameterisation, further experimental work would benefit the robustness of the model greatly. The parameter values for glucuronidation and oxidation pathways are obtainable from the literature, whilst that of sulphation is less well characterised. Though the analysis in this paper suggests that acetaminophen metabolism via the sulphation pathway is secondary in humans, it appears to be important in rats (Reith et al., 2009), which are much more resistant to APAP at human toxic levels. Consequently the critical concentration expression, Eq. (18), will be completely different for rats; we expect that the model is suitably generic to describe acetaminophen metabolism in other species with little modification. However, to fully understand the contrast between rat and human models, for example, more data on metabolism via. sulphation and subsequent model reparamaterisation for the different species is essential.
Through numerical, sensitivity and asymptotic analysis we have improved our understanding of how the different pathways behave. We have highlighted key parameters that our system is sensitive to and also found how the pathways interact with each other, and how this affects the production of protein adducts and the potential for toxicity. This work will provide a foundation on which to build by working directly with scientific researchers and provides us with new areas to research and expand upon using the existing model. This research is part of a larger project funded by the National Centre for the Replacement, Refinement and Reduction of Animals in Research (NC3R's) which aims to improve in vitro testing and reduce the animal testing in science (Krewski et al., 2010). From the initial results and insights, this model is an encouraging first step towards the long-term goal of combining modelling and experimental approaches to mitigate the use of animal testing in toxicological studies, for example, testing hypotheses which would normally be tested in animal models. Its simplicity and analytical tractability means that we can draw conclusions on key parameters that can then be found from in vitro data.