A Combined In Vitro/In Silico Approach to Identifying Off-Target Receptor Toxicity

Summary Many xenobiotics can bind to off-target receptors and cause toxicity via the dysregulation of downstream transcription factors. Identification of subsequent off-target toxicity in these chemicals has often required extensive chemical testing in animal models. An alternative, integrated in vitro/in silico approach for predicting toxic off-target functional responses is presented to refine in vitro receptor identification and reduce the burden on in vivo testing. As part of the methodology, mathematical modeling is used to mechanistically describe processes that regulate transcriptional activity following receptor-ligand binding informed by transcription factor signaling assays. Critical reactions in the signaling cascade are identified to highlight potential perturbation points in the biochemical network that can guide and optimize additional in vitro testing. A physiologically based pharmacokinetic model provides information on the timing and localization of different levels of receptor activation informing whole-body toxic potential resulting from off-target binding.


HIGHLIGHTS
Development of in vitro/in silico framework for identifying off-target toxicity Mathematical modeling of receptor signaling and related transcriptional activity Identification of key events in the signaling pathway Off-target receptor activation in vivo simulated using PBPK modeling

INTRODUCTION
Many drugs are designed to interact specifically with cell surface, cytoplasmic, or nuclear receptors to produce a beneficial therapeutic effect. However, drugs can often bind to and interact with receptors that are not their intended targets, and such ''off-target'' binding may cause what is now often termed a molecular initiating event (MIE), e.g., receptor activation of toxicological relevance that may ultimately lead to an adverse drug reaction (ADR) (Edwards and Aronson, 2000;Guengerich, 2011;Muller and Milton, 2012). In many instances, ADRs can lead to significant morbidity and mortality as well as contribute to high levels of attrition during drug development (Lazarou et al., 1998;Pirmohamed et al., 2004). This can primarily be attributed to an incomplete understanding of the molecular mechanism of action of a given compound and the lack of ability to predict which receptors may be activated unintentionally.
The sole use of in vitro-based experimental strategies in the early stages of drug development and chemical testing is important but can lead to an unreliable and incomplete understanding of reactions (Coleman, 2011). Therefore, often considerable numbers of animals are used to screen out chemicals that may cause off-target toxicity, with figures for the UK reporting that 306,000 in vivo toxicology safety procedures were performed in 2014 (Home Office, 2015). In addition, the chemical industry used almost 345,000 animals in the EU for toxicological or other safety evaluations (European Commission, 2013), and in the United States 3-6 million fish are used annually for whole effluent toxicity testing (Scholz et al., 2013). Furthermore, pharmacokinetics and pharmacodynamics are significantly different between animal models and humans, diminishing their effectiveness in detecting toxicity through pre-clinical studies (Lauschke et al., 2016). There is therefore a clear need to develop scientific approaches to identify toxicologically relevant off-target receptor binding to reduce the burden of animal use in toxicity testing. The development of a more ethical, non-animal toolkit for initial chemical toxicological assessment using an integrated human-based in vitro/in silico system would enhance current strategies and may even expedite the drug development pipeline.
In intracellular signaling, ligand/receptor interactions lead to the activation of a distinct set of transcription factors, the effects of which tend to be tissue specific. Several companies now offer transcription factor activation profiling platforms, and so it is possible to identify and catalog the transcription factor activation profiles of toxicologically relevant receptors upon binding of their known ligands/drugs. It is assumed that transcription factor profiles generated from off-target receptor activation of any given drug can be matched against known ligand/receptor transcription profiles to predict which specific receptor (or class of receptors) has been activated in the initial off-target MIE. However, when testing off-target profiles of new compounds, the resulting transcription profile may not precisely match that of a known receptor (e.g., partial agonism or the binding of multiple receptors), and therefore a method of refinement is required to narrow the subset of off-target receptors. Our approach aims to refine the in vitro receptor identification process for off-target receptors by using information about the changes in receptor-mediated transcription factor activity following the introduction of a given compound and integrating this information with predictive in silico models and analysis. This approach allows for the identification of relevant perturbations in the transcription factor signaling pathway that signify the binding of a receptor or smaller range of receptors as well as other points of interest in the transcription factor signaling network that can contribute toward and guide subsequent off-target receptor identification.
Translating the wealth of knowledge on network interactions of cellular components to dynamic models is generally limited by the amount of available quantitative information to accompany these relationships, such as molecular amounts and reaction rates. However, qualitative dynamic network modeling can be used to compare with routinely generated semi-quantitative experimental time course data, where perturbations can provide valuable information about the system. In silico modeling of this type then provides a platform for the refinement of more quantitative (parameter-based) modeling (Fisher et al., 2013). In such a scenario, the network modeling method of Petri nets provides an effective tool, particularly in the complex, stochastic framework of molecular biological pathways (Chaouiya, 2007;Heiner et al., 2008;Heidary et al., 2015). Petri nets are often used to model multiple species and reactions without defining large quantities of unknown parameters, as modeling emphasis is on network topology and relative amounts of species rather than on specific reaction rates. This emphasis on network structure can then be translated to methods such as flux balance analysis and metabolic control analysis (MCA) without knowledge of rate constants, as was shown for the switching of the metabolic pathway in E. coli (Edwards et al., 2001;Kitano, 2002).
The identification of off-target receptor binding alone for a given compound is insufficient to predict significant off-target toxicity, and so we aim to provide additional information to support and refine the subsequent evaluation of toxic potential. This is achieved by translating knowledge of receptor binding properties and relative distribution of the receptor throughout the body to a whole-body response to the xenobiotic. This approach utilizes a physiologically based pharmacokinetic (PBPK) model adapted specifically for describing receptor activation throughout the body following compound exposure. A PBPK model is a mechanistic, multi-compartment mathematical model that describes the time course dynamics and overall kinetics of an administered drug dose throughout the organism of interest. PBPK models integrate the physicochemical properties of the substance with the specific physiology of the organism such that the evolution of the ADME (absorption, distribution, metabolism, and excretion) processes can be simulated in silico. Drug/substance properties include tissue affinity, membrane permeability, enzymatic stability, etc., whereas the organism/system components include properties such as organ mass/volume and blood flow (Rowland et al., 2011). PBPK modeling is used in this work to couple the pharmacokinetics of a drug to dose-response parameters with the associated off-target receptor in different tissues to generate spatiotemporal dynamics of the off-target receptor activation.

Development of the Signaling Pathway Model
As proof of concept, an in silico model of the histamine H1 receptor signaling pathway was formulated. This pathway was chosen owing to the well-understood intracellular signaling interactions involved upon receptor stimulation and the existence of a known off-target partial agonist, lisuride (Bakker et al., 2004). The H1 receptor is a G-protein-coupled receptor that, upon activation, leads to dissociation of Ga q/11 and the Gbg complex. Ga q/11 activates phospholipase Cb (PLCb), leading to hydrolysis of phosphatidylinositol 4,5-biphosphate (PIP 2 ) and the formation of inositol triphosphate (IP 3 ) and diacylglycerol (DAG) (Bakker et al., 2001;Sandal et al., 2013). IP 3 mediates a transient intracellular calcium release from the ER (Shah et al., 2015) that eventually mediates the activation of nuclear factor of activated T-cells (NFAT) (Macian, 2005), cAMP response element-binding protein (CREB) (Johannessen and Moens, 2007), and myocyte enhancer factor-2 (Mef2) transcription factors (Lu et al., 2000). DAG simultaneously activates protein kinase C (PKC), and this phosphorylates IkB kinase (IKK), ultimately leading to nuclear factor kappalight-chain-enhancer of activated B cells (NF-kB) transcription factor activation (La Porta and Comolli, 1997). The Gbg complex also plays a role in histamine signal transduction, regulating many effectors, including adenylate cyclase (AC) (Maruko et al., 2005) and phosphoinositide 3 kinase (PI3K) (Gautam et al., 1998). AC mediates the subsequent activation of protein kinase A via cyclic AMP (cAMP) leading to CREB phosphorylation and transcription factor activation (Mosenden and Taské n, 2011). PI3K mediates the activation of Akt, NF-kB, and activating transcription factor 2 (ATF2) (Bence et al., 1997;Breitwieser et al., 2007).
To provide semi-quantitative information for the relative transcription factor dynamics as described earlier, we assayed pathway perturbations using a luciferase reporter-based transcription factor array to calibrate the fold increase expected of key signaling outputs upon stimulation with an agonist. These transcription factors were identified as NFAT, NF-kB, CREB, Mef2, and ATF2. Incubation of H1 receptor-expressing HeLa cells with histamine showed considerable activation of these transcription factors (Table 1).
A stochastic Petri net model of the histamine H1 receptor signaling pathway was formulated based on existing knowledge of the pathway and network interactions with the five critical transcription factors determined to be activated following ligand binding. The pathway in this proof of concept provides an illustrative example of what should ultimately form part of a larger cell signaling model that incorporates the complexity of the known toxicological receptors and associated transcription factors in the proposed methodology. The H1 Petri net includes the key dynamic molecular species and appropriate network interactions that are activated during ligand-binding-induced signaling. This pathway is depicted using the modified Edinburgh pathway notation (mEPN) format (Freeman et al., 2010) in Figure 1 and directly corresponds to the layout of the Petri net. All rates are equal such that all stochastic transitions are equally likely to fire but are effectively modulated by the concentration of upstream reactants in a mass action process. Time is interpreted qualitatively reflecting the relative order of events. Varying quantities in the mathematical model, such as the amount of ligand introduced (''dose'') and the total amounts of system species (i.e., moieties of active and inactive states for each protein), modulate the scale of transcriptional activity regulation, and as such, these values were optimized to correlate with the experimental signaling assays. This optimization was carried out by assuming a large-scale continuum approximation of the Petri net to a system of ordinary differential equations (ODEs) and fitting to the corresponding transcription factor output data ( Figure 2). It should be noted that the optimal parameter set is non-identifiable for such a large system with relatively few data points to fit. However, this issue was the precise motivation for the combined Petri net/MCA approach, which is well suited to understanding the relative impact of small perturbations on the transcription factors of interest and prioritize network connectivity information in favor of accurate predictions of parameters and dynamics (Koch et al., 2010). Corresponding pathway reactions, moieties, and ODEs can be found in the supplementary material. In addition to providing static information on the network interactions of the signaling pathway and relative changes in steady state activity following receptor activation, Petri nets can also be used to simulate transient temporal dynamics, providing further dynamic information on the relative order and scale of transcriptional regulation ( Figure 3) following a receptor-ligand binding event. However, it is clear that more data would be required for one to relate this dynamic output to the biological context and validate any potential predictions about transient dynamics.

Analysis of Network Perturbations to Identify Off-Target Responses
The identification of significant pathway reactions upstream of transcription was achieved using metabolic control analysis (MCA), which is a mathematical technique that tests the sensitivity of a given variable to network perturbations (Kacser and Burns, 1973;Heinrich and Rapoport, 1974 concentration control coefficients provide the ratio between a relative measure of change in the steady state of a system variable as affected by perturbations in network reaction rates. In our illustrative H1 example model, MCA coefficients were calculated for each transcription factor that was experimentally determined to show significant change in activity following binding of the H1 receptor ( Figure 4). The rows of the heatmap in Figure 4 correspond to the numbered reactions as indicated in the supplementary material. MCA points not only to the direct regulation of gene transcription as critical to H1-associated The Petri net describes the key relationships between components of the signaling pathway system culminating in the regulation of downstream transcription factor expression stimulated by the binding of a ligand to the histamine H1 receptor.
transcriptional activity (white patches in Figure 4), but also to other reactions within the cascade, upstream of the transcription factors and downstream of the target receptor. For example, in this system the transcriptional activity of Mef2 is sensitive to relatively distant biochemical reactions, such as the rate of calcium release from the ER (24% of maximum sensitivity provided by perturbation of Mef2 transcription rate). Also, the model suggests that the transcriptional activity of ATF2 is more sensitive to perturbations in PIP2 synthesis than to the regulation of the BTK:PIP3 complex that directly activates ATF2 by phosphorylation.
The identification of these sensitive perturbation points within the signaling pathway model provide information beyond the transcription factor activity measurements found experimentally, which allows for more optimized, directed experimental designs for receptor identification, if initial screening fails to identify the off-target receptor. For example, for a given compound that was shown to regulate Mef2 transcriptional activity but did not interact with the H1 receptor, this model would inform a proposal to screen for receptors that are known to interact with biochemical reactions identified as being sensitive, such as calcium release, during MCA.

Translation to Tissue Scales Using a PBPK Model
Following an in silico identification of an off-target receptor, extrapolation to the study of potential in vivo toxicity can be performed using a PBPK model. For our illustrative example, receptor binding properties are provided by EC 50 dose-response curves for the off-target H1 agonist, lisuride ( Figure 5A), and measurements of the corresponding binding affinity, K d (Bakker et al., 2004). The dose-response curves were estimated by fitting the following equation to the dose-response data: for ligand concentration L. The optimized parameter values are given in Table 2. To provide tissuespecific responses we also used western blot measurements of relative H1 receptor expression in different tissues (Figures 5B and 5C) and calculated modified tissue-specific EC 50 values using where i denotes the i th tissue, K d is the dissociation equilibrium constant for lisuride, and R i is a measure of receptor abundancy in tissue i (see Table 3). For simplicity, this model assumes that the same amount of receptor binding is required to achieve 50% response in each tissue in the absence of any other information, particularly as the response measured is proximal to receptor binding attenuating any potential amplification effects arising from potential signaling cascades in different tissues (Kenakin, 2009). For further information regarding this derivation see the Supplemental Information.
To simulate the pharmacokinetics of lisuride throughout the body, physicochemical properties of the compound were required, which were obtained from previously published measurements. These properties include lipophilicity, whether the drug is neutral/acid/base, solubility (obtained from the DrugBank database [Wishart et al., 2006]), molecular weight (O'Neil, 2013), acid dissociation constant (Meloun et al., 2005), and effective permeability (Winiwarter et al., 1998). The time course dynamics simulated by the PBPK model for drug concentration in each tissue compartment of the body were then coupled to receptor binding properties and relative receptor expression in tissues to provide a predictive temporal response throughout the body. This response can be produced for any dosage regimen and various methods of administration, such as intravenous, oral, and inhalation. The PBPK model was based on the form derived by Peters (2008). The model was optimized for lisuride physicochemical and binding properties and the H1 receptor distribution throughout the different tissues. Example lisuride response kinetics following both intravenous (IV) and oral administrations can be found in Figure 6. The IV dose of 25 mg/mL used in Figure 6 was the same as that used in a previous pharmacokinetic study for relevance (Krause et al., 1991). These experimental data were also the IV data used to optimize the PBPK model to recapitulate the lisuride dynamics in the venous blood compartment and also simulate corresponding oral profiles as per the methodology described by Peters (2008). The oral dose of 0.1 mg chosen for the PBPK model was deemed relevant by matching previous pharmacological studies ( Koizumi et al., 1985;Al-Sereiti and Turner, 1989). The dynamic response of the H1 receptor is visualized over time as a solution to Equation (1 with tissue-specific EC 50 values for the pharmacokinetics of lisuride (L) in different parts of the body. Both IV and oral administration simulations are plotted to also highlight the impact of delivery route. This is particularly pertinent in this case where we are studying a receptor that has a relatively high concentration in the gastrointestinal tract. IV administration results in relatively high receptor stimulation in the liver, brain, small intestine, and colon at earlier times, whereas oral administration results in a more gradual accumulation in these tissues and the receptors in the colon are stimulated at a near-maximal level for a relatively long time after oral ingestion. These simulations allow us to compare how the off-target response varies throughout the body over time depending on the pharmacokinetics of the drug coupled with physiologically relevant receptor availability and receptor binding information. Such information is potentially useful to determine whether or not an identified off-target agonist is likely to elicit an off-target receptor response in an area of high target density based on its physicochemical properties.

DISCUSSION
ADRs are a major cause of patient morbidity, mortality, and drug attrition during development (Pirmohamed et al., 2004). This can be attributed to a poor understanding of the mechanisms underlying the toxic response and also to a lack of current tools for the prediction of a toxic outcome. Animal models have a limited scope, and data obtained using such models may not be ideal for ascertaining toxicity seen in humans. As such, computational systems biology models can be essential tools to improve chemical reaction predictivity (Krewski et al., 2010). In this study, we describe a new in silico modeling method that can be used to enhance the current knowledge of pathway perturbations to provide a new toxicity-testing paradigm based on human biology. In this method, chemical-mediated activation of transcription factors and intracellular signaling pathway molecules were used as readouts to inform and drive a pathway-based in silico approach to identify possible upstream receptor(s) engaged by such chemicals. In vitro data were then used to inform a PBPK in silico modeling platform to understand and rank the risk of toxicity at tissue, organ, and whole-body levels over time. Key to this integrative approach was the coupling of in vitro experimental techniques and advanced in silico modeling to create a unique resource that, with further  development and parameterization, could be used to predict the off-target toxicity of compounds that can then inform and direct more focused in vivo experimentation.
Mathematical modeling was used to mechanistically describe the processes that lead to regulation of transcriptional activity following the binding of ligand to receptor. This was achieved by designing a signaling pathway model that represented all the relevant processes and biochemical reactions downstream of ligand binding, culminating in the regulation of transcription. We have established a novel in vitro/in silico approach using data from assays measuring transcription factor activation and chemically induced perturbations of intracellular signaling pathways to inform in silico pathway modeling. This unbiased pathway-led approach uses computational simulations to identify causality between receptor activation and pathway perturbations to aid identification of the upstream receptor/s engaged by the initial MIE. As proof of concept, an in silico Petri net model of the histamine H1 receptor-signaling pathway was formulated with the off-target compound, lisuride. The output of this system provides semi-quantitative temporal dynamics for the entire pathway that can be used to investigate system perturbations, simulate experiments, and provide structural pathway predictions. In vitro reporter assay data were then used to parameterize and validate the model, and the identification of critical candidate perturbation points was achieved using MCA. Signaling pathway models can be purposely used in this methodology to provide a library of MCA coefficients for a range of transcription factors associated with receptor binding and toxicity and guide further experimentation. In the example shown, calcium release from the ER and PIP2 synthesis are highlighted as important upstream events for the transcriptional activity of Mef2 and ATF2. If a new compound is shown to induce the activity of these transcription factors but the receptor responsible is not identified via screening, for instance, further testing could be guided toward targets that modulate these upstream processes. This illustrates the feasibility of this approach in directing further experimentation toward relevant pathway mechanisms or receptor clusters during the process of receptor identification via focused in vitro assay testing.
In vitro to in vivo extrapolations of whole-body consequences of receptor binding was explored using PBPK modeling. The structure of PBPK models typically revolves around the anatomical structure of the organism, with different organs and tissues of varying perfusion rates being separated into distinct compartments. These compartments are then coupled through the circulation, whose arterial and venous flow is described to connect the organs in a physiological way. Entrance points (e.g., absorption) of the model depend on the drug administration method (e.g., inhalation, ingestion, injection), whereas exit points (e.g., excretion) are generally described via the kidneys and intestine. The flow kinetics of the model determine distribution, whereas metabolism occurs in the liver and intestine. The inherent physiological basis distinguishes true PBPK models from their PK model counterparts that usually simplify the physiology to fewer hypothetical compartments of different flow rates, driven by the data/process of interest, such that they are often more tractable analytically. In contrast, PBPK models are generally more complex but are designed to have a better global representation such that valid extrapolations can be made and disparate experimental data can be integrated during model parameterization. In this way, PBPK models are less reliant on data fitting to obtain appropriate values for equation parameters and essentially the same model (with appropriate modifications) can be suitably applied in many different pharmacological scenarios for quantitative risk assessment and therapy optimization.
PBPK model simulations are increasingly being used in pharmacology, in both academia and industry, to provide important predictions of the pharmacokinetic properties and toxic potential of new drugs at an early stage in drug development (Zhao et al., 2011;Jones and Rowland-Yeo, 2013;Tsamandouras et al., 2015). This type of in silico testing can offer a quicker, cheaper, and more ethical alternative method when compared with traditional in vivo experiments performed. Ideally, both experimental and computational methods are used harmoniously to provide a cycle of information and enhanced knowledge iteration as the accuracy of PBPK models inevitably rely on quality experimental data to calibrate rates within the differential equations. In the method reported here, physicochemical properties of the chemical are combined with tissue-specific receptor expression and EC 50 data to predict time course dynamics of the chemical concentrations in each tissue, as well as tissue-level receptor-activation responses to that chemical. These predictions can be produced for any dosage regimen and various methods of administration. In the example study of the off-target partial agonist of the histamine H1 receptor, lisuride, the combination of lisuride pharmacokinetics and relative H1 receptor distribution throughout the body allowed us Receptor activation of the H1 histamine receptor was studied with known agonist (histamine) and off-target agonist (lisuride). Using these assays, each parameter was calculated using GraphPad Prism.
to predict that the dose response would be most significant in the brain, liver, and gastrointestinal system. In this case example, these results are supported by prior knowledge of the compound and receptor, although the modeling was done agnostic of such prior in vivo findings. In particular, receptor response localized to the brain is somewhat expected since lisuride is primarily a psychotherapeutic drug, affecting dopamine and serotonin regulation (Marona-Lewicka et al., 2002). Lisuride is primarily metabolized in the liver, where there is a relatively high expression of histamine receptors. There is also high receptor expression in the gastrointestinal tract owing to the role of histamine in intestinal secretion and motility (Leurs et al., 1995;Sander et al., 2006). Furthermore, lisuride administration in patients with Parkinson disease has been associated with gastrointestinal side effects (Ebadi and Pfeiffer, 2004). Although relative response rates have been quantified by the model in different parts of the body at different times, to translate what such a response directly represents in the context of toxicity and clinical relevance is very complicated, and restricted in this methodology, establishing a challenge beyond the scope of this paper. However, these PBPK-based extrapolations do allow us to generate predictive data relevant to risk assessment and further translation to toxicity at the organ and whole-body levels for off-target receptor perturbations. The output provided by this method is intended to identify toxic potential and guide subsequent in vitro and in vivo experimentation to organs of interest/importance.
The operating parameters of the approach are circumscribed by the extent of current knowledge regarding receptors and their function. This represents a potential limitation of the strategy, although the mathematically driven signaling pathway model has the potential to identify novel, uncharacterized receptor targets. The challenge of identifying sensitive perturbation points within large-scale networks of receptor signaling pathways required that a semi-quantitative network-based approach be used. This inevitably limits the amount of predictive, dynamic information that can be extrapolated, and caution must be exercised such that the utility of mathematical models is preserved by acknowledging the relevant application that stimulated its design. The approach is experimental (with elements of modeling and extrapolation to assess and rank toxicological risk) and does not incorporate prediction of receptor binding based on chemical or receptor structures. The strength of the methodology is predicated on currently available, validated experimental methods, as it does not require the development of new, untested technologies and relies on sound criteria-based selection of receptors and quantifying receptor function and binding using established experimental techniques. Future work requires the development of multiple pathway models based on training chemical data as well as the integration of pathways, which should be optimized and validated with non-training data. Furthermore, the current PBPK framework can be extended to ensure improved predictive potential by incorporating mechanistic tissue models, catering for a wider range of chemicals and capturing population-level responses. More work is also needed to translate tissue-level receptor activation responses to measures of toxicity, such as relevant biomarkers. Carefully calculated person-to-person variation and covariances within organism-related parameters would also allow for the prediction of a population response whereby different individuals within a sample population may exhibit different levels of exposure and therefore associated toxicity from the same dosage levels.  approach of this study has shown how the multidisciplinary, iterative process of systems biology can be applied to direct experiments, optimize the utility of generated data, and challenge and refine theoretical modeling to improve methods for detecting and predicting toxicity caused by compounds that bind to offtarget receptors.

METHODS
All methods can be found in the accompanying Transparent Methods supplemental file.

Cell culture
HeLa cells were cultured in Dulbecco's Modified Eagle Medium supplemented with 10% foetal bovine serum, 2 mM glutamine and incubated at 37 °C in a 5% CO 2 humidified atmosphere.

Quantification of transcription factor activation
Transcription factor activation was measured using Cignal Reporter dual-luciferase assay (Qiagen) as per manufacturer's instructions, to determine intracellular signal transduction perturbation and transcription factor activation data for in silico modelling. Briefly, HeLa cells were plated in 96-well plates and transfected with inducible transcription factor responsive constructs composed of specific transcription factor response elements linked to firefly luciferase along with a constitutively expressing Renilla construct as an internal control. Transfection was performed with Attractene transfection reagent (Qiagen) as per manufacturer's instructions. Transfected cells were stimulated with histamine (Sigma) for 6 or 24 h and activation of specific transcription factors quantified by luminometer using the Dual-Luciferase Reporter Assay System (Promega).

Receptor activation
Activation of histamine receptor H1 with its known agonist, histamine and the off-target partial agonist, lisuride was determined using PathHunter® Chinese Hamster Overy (CHO)-K1 reporter cells that over-express ProLink™-tagged human histamine H1 receptor and the Enzyme Acceptor (EA)-tagged β-Arrestin (DiscoverX) as per manufacturer's instructions. In these cells activation of histamine receptor H1 by a ligand forces H1 receptor coupling with EA-tagged β-Arrestin leading to the formation of a functional enzyme that is able to generate a chemiluminescent signal upon substrate cleavage. These assays were used to calculate EC50 values for both histamine and lisuride.

Immunoblotting analysis
Mouse organs from wild type C57 mice (Charles River) were snap frozen in liquid nitrogen and whole tissue lysates homogenised in lysis buffer (250 mM sucrose, 50 mM Tris.HCl,1 mM MgCl 2 , 1% Triton X-100, 1 mM PMSF). Lysates were also generated from HeLa and CHO cells over-expressing histamine H1 receptor. The lysates were centrifuged at 5000 g for 15 minutes at 4 °C and 20 µg of total protein separated using SDS-PAGE and analysed using standard immunoblotting techniques. Membranes were blocked overnight in 5% milk powder in blocking buffer (0.2% Triton X-100 in PBS) at 4 °C and then incubated with rabbit polyclonal anti-H 1 R antiserum (Abcam; 1:1000) followed by horseradish peroxidaseconjugated goat anti-rabbit antiserum (Jackson Labs; 1:2500). Immuno-reactive proteins were visualized using enhanced chemiluminescent substrate detection (ThermoFisher). Densitometry analysis was carried out using Image J and fold-changes in expression relative to HeLa cells calculated.

Petri Nets
The stochastic Petri net model of the histamine H1 receptor signalling pathway was formulated using Snoopy software (Heiner et al., 2012). In the context of cellular signalling pathways, a Petri net is a bipartite graphical representation of a biochemical network with "tokens"/dots-or-numbers (representing number of molecules or concentration) in "places"/circles (network species, proteins etc.). A reaction is represented by the "firing" of a "transition"/square where tokens are moved from one place to a downstream place as indicated by connecting "arcs"/arrows. A transition may fire only if all upstream places contain sufficient tokens, i.e. all reactants must be present for a specific reaction to occur. Petri nets use mass action and Gillespie algorithm principles such that transitions with a higher number of upstream tokens fire with a higher probability.

Parameter Optimisation and Metabolic Control Analysis
The signalling pathway model large-scale continuum approximation to an ODE system was solved using a Runge-Kutta 4/5 th order method with Matlab R2015b software. Pathway moieties (conserved quantities defined in supplementary material S1.2) were optimised such that model signalling output was consistent with experimental data. Optimisation was carried out by minimising the error between fold-increase in transcription factor levels in the model and the data. Metabolic scaled concentration control coefficients were calculated using Copasi 4.15 software (Hoops et al., 2006). The Peters (2008) PBPK model code was written and solved with Matlab R2015b software.
Table S1, related to Figure 1: Histamine H1 receptor signalling pathway reactions. Each reaction in the Petri net of Figure 1 is described explicitly with associated description and mass action kinetics.
Pathway moieties were optimised to fit the steady state fold-changes of the transcription factors ( Figure 2 of the main manuscript). Due to the nature of the approach taken (maximising network/pathway connectivity information rather than focusing on dynamics with a minimal model), there is an inevitable issue with model parameterisation. We have made efforts to find the global optimum for the parameter set through Latin hypercube sampling, but it is clear that other parameterisations could give fits that would at least look as good (by eye) and fit well within any expected variation from the biology/experimental error.
A profile likelihood estimation method was employed to determine parameter identifiability analysis and, as expected, the results indicate that the parameters are not uniquely identifiable for such a problem since the number of free parameters far outnumber the number of data. The analysis was performed using Data2Dynamics software (Raue et al., 2013, Raue et al., 2015 and the profile likelihood method to analyses parameter identifiability, as developed by Prof Jens Timmer's group in Freiburg (Raue et al., 2009, Maiwald et al., 2016. The results of this identifiability analysis are illustrated in Figure S2. Additionally, local sensitivity analysis indicates that, in line with the complementary MCA method (as described in the main manuscript), it is the nodes that are relatively close to the transcription factor signal that dominate the subsequent expression following receptor activation. These fairly intuitive results suggest that errors arising due to parameter changes throughout the system can be mitigated for in order to achieve the same fold-changes seen experimentally, provided that (a) other parameters are also changed (in the case of sensitive proximal nodes) or (b) because their values are not having a large effect on the downstream signal (distal nodes). This emphasises the importance of the use of MCA in the proposed methodology as a means of identifying any distal nodes that do affect transcription factor signal, as identified in the paper. Of course, to acquire identifiable parameters, far more data would be needed (than is available for such a complex network) or typically, a minimal ODE model would be constructed that significantly reduces the number of variables and parameters that cannot be identified with given data. This minimal approach would potentially be useful for answering questions about the dynamics of the H1-histamine signalling pathway for a very specific scenario and/or ligand. However, this approach would be limited when trying to identify off-target receptor toxicity. Figure S2, related to Figure 2: Identifiability analysis for the H1 Histamine pathway Petri net. Briefly, following the identification of a globally optimal parameter set each model parameter is perturbed within a prescribed range, one parameter at a time. Starting from the optimal value, a single parameter is perturbed, set to be fixed at the new and non-optimal value, and the model is then re-fit, i.e. by optimising all other parameters. Then corresponding likelihood values are plotted for this new fit. Analysis was conducted using the profile likelihood method in Data2Dynamics software, solved in MATLAB 2017b.