Email updates

Keep up to date with the latest news and content from BMC Systems Biology and BioMed Central.

Open Access Methodology article

In silico labeling reveals the time-dependent label half-life and transit-time in dynamical systems

Thomas Maiwald12*, Julie Blumberg34, Andreas Raue1, Stefan Hengl1, Marcel Schilling5, Sherwin KB Sy6, Verena Becker257, Ursula Klingmüller57 and Jens Timmer11089

Author Affiliations

1 Center for Systems Biology, Freiburg, Germany

2 Department of Systems Biology, Harvard Medical School, Boston, USA

3 Epilepsy Center, University Hospital Freiburg, Germany

4 Department of Neuroscience, Children's Hospital, Harvard Medical School, Boston, USA

5 Division Systems Biology of Signal Transduction, DKFZ-ZMBH Alliance, German Cancer Research Center, Heidelberg, Germany

6 College of Pharmacy, University of Florida, Gainesville, USA

7 Bioquant, Heidelberg University, Germany

8 Freiburg Institute for Advanced Studies, University of Freiburg, Germany

9 BIOSS Centre for Biological Signalling Studies, University of Freiburg, Germany

10 Department of Clinical and Experimental Medicine, Linköping University, Sweden

For all author emails, please log on.

BMC Systems Biology 2012, 6:13  doi:10.1186/1752-0509-6-13

The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1752-0509/6/13


Received:10 May 2011
Accepted:27 February 2012
Published:27 February 2012

© 2012 Maiwald et al.; licensee BioMed Central Ltd.

This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

Mathematical models of dynamical systems facilitate the computation of characteristic properties that are not accessible experimentally. In cell biology, two main properties of interest are (1) the time-period a protein is accessible to other molecules in a certain state - its half-life - and (2) the time it spends when passing through a subsystem - its transit-time. We discuss two approaches to quantify the half-life, present the novel method of in silico labeling, and introduce the label half-life and label transit-time. The developed method has been motivated by laboratory tracer experiments. To investigate the kinetic properties and behavior of a substance of interest, we computationally label this species in order to track it throughout its life cycle. The corresponding mathematical model is extended by an additional set of reactions for the labeled species, avoiding any double-counting within closed circuits, correcting for the influences of upstream fluxes, and taking into account combinatorial multiplicity for complexes or reactions with several reactants or products. A profile likelihood approach is used to estimate confidence intervals on the label half-life and transit-time.

Results

Application to the JAK-STAT signaling pathway in Epo-stimulated BaF3-EpoR cells enabled the calculation of the time-dependent label half-life and transit-time of STAT species. The results were robust against parameter uncertainties.

Conclusions

Our approach renders possible the estimation of species and label half-lives and transit-times. It is applicable to large non-linear systems and an implementation is provided within the PottersWheel modeling framework (http://www.potterswheel.de webcite).

Background

Motivation

An increasing number of biological phenomena are described by mathematical models, specifically on the basis of biochemical reaction networks [1,2]. The dynamic properties of these networks are given by their model structure, kinetic parameters, initial values of the involved species, and externally specified input functions. The interpretation of an isolated element of the network, e.g. a certain rate constant, has only a limited meaning, because its effect can only be understood when taking the whole network context into account. We therefore seek to introduce two dynamical characteristics which have a physiological meaning, are intuitive to understand, and capture the system kinetics on a higher level of abstraction. The first characteristic, the label half-life, applies the half-life concept not to a species, but to a virtual label attached to the species. The second one, the label transit-time, is the time-period it takes for a fraction of labeled entities to pass through a subsystem of the network. Both quantities are calculated using a novel approach called in silico labeling, which is also introduced in the present work.

In Silico Labeling and Species vs. Label Half-Life

In a laboratory tracer experiment, a substance is marked to better understand the kinetic properties of the dynamical system [3]. Different tracer substances have been used, e.g. radioactive iodine-125 [4,5] or green fluorescent protein-tagged proteins in combination with fluorescence recovery after photobleaching (FRAP) [6]. A good tracer does not hamper the flux of the substance, therefore one can assume that the flux of the tracer within a certain reaction is proportional to the flux of the original species. This is the key property of the in silico labeling approach, where an additional set of reactions is added to an existing mathematical model describing the kinetic behavior of a tracer, called the label. In contrast to real tracer experiments, the in silico method offers the opportunity to define dead-ends, avoid double-counting of cycling label, and to restrict the label to a sub-network of reactions. This allows asking specific questions about the original system, like how long it takes for 50% of the molecules of a substance to travel along a certain path, while in reality an alternative path may exist. In addition, predominant paths can be identified in deterministic models as has been done previously for stochastic systems [7].

Mathematically, the half-life T1/2 of a species is defined as the time-period until it reaches half of its initial amount assuming no influx. For clarity, we denote this time-period as the species half-life (SHL). In non-isolated and non-linear processes, this time-period differs from the amount of time required for 50% of initially existing molecules to be processed. For this, we introduce the label half-life (LHL), defined as the half-life of the label of a species. Equalities and differences between the species and label half-life are displayed in Figure 1 and proven in the methods section.

thumbnailFigure 1. Species vs. label half-life. Panel A: The species half-life of a substrate S in the reaction S P is plotted for different reaction types (solid lines). Except for processes of order 1, the half-life is time-dependent. Since the substrate is not produced in further reactions, the label half-life (dashes) equals the species half-life. Panel B: The substrate S participates in a production (A S) and a processing (S P ) reaction. Now, the species half-life differs except for a linear processing from the label half-life, because the label flux is proportional to the total flux of each reaction and is therefore affected by concentration changes through influx of S. Both panels: The species half-life has been determined analytically and numerically according to the methods section. Matlab scripts to reproduce the plots are available in the additional file 1.

While for simple systems the species half-life can be determined analytically, the symbolic integration of a Michaelis-Menten kinetics leads to advanced mathematical calculations including the Lambert W function [8]. We therefore also provide an automatic and generally applicable numerical method to determine the species half-life.

Label Transit-Time

Transit-times are discussed in a variety of fields and they are, for example, used to quantify how quickly food moves through the gastrointestinal tract [9]. When describing the dynamics of Markovian particles, the mean transit-time denotes the time spent on average in a subsystem [10], while the mean sojourn-time also takes into account the probability that the subsystem is entered at all [11]. In pharmacokinetics, the so-called mean residence time values [12] are estimated based on empirical data assuming linear kinetics [13]. Apart from linearity, no influx for the species of interest is permitted. Eventually, the estimation is only applicable to observable species. The computation of the mean residence time is accomplished by the ratio of the area under the first moment curve (AUMC) to the area under the curve (AUC) of the concentration-time profile of a drug [14].

We here introduce the label transit-time (LTT) from a source to a target pool in a chemical reaction network as the time-period after which 50% of all entities residing in the source pool at t = 0 have reached the target pool at least once. The exact path from source to target pool is not important in the unconditioned case. The LTT information could be valuable to estimate the time for a drug or an enzyme to reach its site of action.

Extended Reaction Network

To determine the label half-life, it is important to distinguish entities residing in the source pool at t = 0 from other entities entering the source pool at later time-points. When calculating transit-times, this discrimination has to be applied to all pools and fluxes between source and target. To achieve this aim, the species of interest is computationally labeled and subsequently tracked throughout the dynamical model. The labeling is realized by an additional set of reactions describing the kinetic behavior of the labeled species, depending on the kind of time characteristic LHL or LTT, the source species, and potentially a target species.

In case of label half-life calculations, it is sufficient to create labeled reactions for all reactions in which the source species is a reactant. In fact, labeled reactions are prohibited if the source species is a product; this is to avoid double-counting the labeled species. In the case of transit-time calculations, for all original reactions in which labeled species are involved, a new labeled reaction is added. In all labeled reactions with the target species being the product, the label is removed and accumulated in an artificial pool which is used to determine when 50% of the existing label has reached the target.

The label stays virtually attached to a species throughout all modifications of the species, such as phosphorylation or relocalizations, e.g. shuttling into the nucleus. While the suggested approach can be implemented in a straightforward way for monomeric reaction networks with only up to one labeled reactant and product, for the general case where the reactions involve multiple reactants and products or where labeled species may form a polymer, a systematic book-keeping of all possible combinations of labeled and unlabeled species is required.

As motivated by laboratory tracer experiments the fluxes of the additional system are based on the corresponding fluxes in the original one, which is explained in detail in the methods section.

Profile Likelihood-based Confidence Intervals

Recently, we suggested a profile likelihood-based approach to determine the confidence intervals on calibrated parameter values in mechanistic mathematical models [15]. The same reasoning can be applied in order to estimate confidence intervals for the time-dependent label half-life and transit-time characteristics.

Implementation

All concepts have been implemented within the PottersWheel modeling and parameter estimation framework that is available from http://www.potterswheel.de webcite[16] and have recently been applied by the authors to the mathematical models of the erythropoietin and epidermal growth factor receptors [17,18]. The application of the method within the PottersWheel framework is described in additional file 1.

Additional file 1. Application within PottersWheel. This additional file contains MATLAB scripts to run various tasks related to the in silico labeling approach. http://www.biomedcentral.com/imedia/4654854926777309/supp1.pdf. webcite.

Format: PDF Size: 328KB Download file

This file can be viewed with: Adobe Acrobat ReaderOpen Data

In the next section, the proposed labeling method is illustrated for the JAK-STAT signal transduction pathway and afterwards described in detail. After proving the equality of species and label half-life for isolated or linear processes, a fitted model of the JAK-STAT pathway is used to determine the label half-life of unphosphorylated STAT and its label transit-time when cycling through the nucleus of a cell.

Methods

Illustration of the method

Figure 2 illustrates the in silico labeling approach for the JAK-STAT signal transduction pathway, where STAT molecules cycle between cytoplasm and nucleus. First, cytoplasmic STAT molecules (S) are phosphorylated (pS) by an active receptor (pR) and form dimers (pS_pS). The complexes enter the nucleus (npS_npS) where they act as transcription factors, disassociate and are dephosphorylated (nS) again. Finally, they return to the cytoplasm (S) and can be activated again. In order to determine the label half-life of cytoplasmic STAT and the label transit-time for a whole cycle, we set source and target species to unphosphorylated cytoplasmic STAT. At t = 0, all molecules of the source pool are labeled, symbolized by the small red spheres. The label is not removed until the target pool is reached, in this case when a STAT molecule leaves the nucleus. Then, the label is accumulated in an artificial pool of returned label and an unlabeled STAT molecule enters the cytoplasm. Over time, the fraction of labeled to free, unlabeled STAT molecules, SL/SF, decreases in the cytoplasm. The total flux <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M1">View MathML</a> of the first reaction, that is the phosphorylation of STAT molecules, can be divided into the flux <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M2">View MathML</a> of labeled and the flux <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M3">View MathML</a> of free STAT molecules. The fraction of the fluxes is set to match the fraction of labeled to free STAT molecules, by the following relationship:

thumbnailFigure 2. In silico labeled JAK-STAT signaling pathway. STAT molecules S (blue) are phosphorylated by an active receptor-kinase complex (pR) and form dimers (pS-pS). These dimers enter the nucleus, dissociate, and are subsequently dephosphorylated. Finally, the single STAT molecules re-enter the cytoplasm, where they can again be phosphorylated and thus continue the nuclear-cytoplasmic shuttling. The labeling approach is visualized by red spheres attached to the STAT molecules. At t = 0, all cytoplasmic STAT molecules are labeled. After the nuclear export, the label is removed from the molecule and enters the artificial pool of returned label. Consequently, an increasing fraction of cytoplasmic STAT molecules are not labeled which has to be considered in the calculation of fluxes for free and labeled entities. To determine the time-dependent label half-life and transit-time values, the labeling procedure is repeated for a series of time-points.

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M4">View MathML</a>

(1)

The label half-life of STAT at time-point t is given by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M5">View MathML</a>

(2)

The label transit-time from STAT to STAT at time-point t can be derived from the time-profile of the returned label RL:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M6">View MathML</a>

(3)

This procedure is repeated for a series of time-points t in order to determine LHL(t) and LTT(t) for all time-point of interest.

Terminology

In the following we assume that the biological system is mathematically described by a set of reactions rj, 1 ≤ j n, corresponding to a set of coupled differential equations. The concentration change of each entity xi, 1 ≤ i ≤ m, is the sum over all fluxes of reactions where it appears as a product minus the sum over all fluxes of reactions where it appears as a reactant, mathematically [19]

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M7">View MathML</a>

(4)

Here, vj describes the flux of reaction j, aij ≥ 0 the stoichiometry of xi as a product in reaction j and bij ≥ 0 the stoichiometry of xi as a reactant in reaction j. We use the same symbol for an entity and its concentration, [xi] ≐ xi. The time-profile of each species can then be calculated for given initial values <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M8">View MathML</a> and potentially driving input functions uk(t). The flux vj of reaction j may be a non-linear function of one or more species concentrations xi and externally defined uk. To improve readability, we omit explicitly denoting the time-dependency, i.e. xi(t) is rather written as xi.

Analytical and numerical half-life calculation

The half-life of a species xi of interest is determined by extending the differential equation network (4) by one equation for an artificial quantity y depending only on the outfluxes of xi,

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M9">View MathML</a>

(5)

with initial value <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M10">View MathML</a> The whole system (4, 5) is solved either analytically or numerically and the species half-life of xi is given by T1/2 for which

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M11">View MathML</a>

(6)

Note that a half-life characterizes the decay of a quantity, independent of any production rates. Therefore, all influx contributions are neglected in equation (5). In general, only linear processes possess a constant half-life. Otherwise, the half-life depends on the initial concentration <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M12">View MathML</a> and is therefore time-dependent. In this case, the above procedure is repeated for a series of different initial time-points t0. In a numerical integration, it is important to limit the maximum integrator step size for an accurate approximation of the y0/2 threshold crossing.

The half-life of a species xi is only partially related to the time it takes for 50% of an experimental tracer to leave the source pool. The two values coincide if xi has either no influx or when the outflux from xi is described by a linear process, which will be proved in the next two subsections. Therefore, we suggest the in silico labeling half-life as a means to determine a time-characteristic which is motivated by laboratory tracer experiment with the additional property to avoid tracer-double counting in kinetic cycles.

In silico labeling half-life for isolated processes

For simplicity, we assume that the species of interest x ∈ {x1, . . . , xm} is consumed only in one reaction. In in silico labeling, the flux of the corresponding label z depends on the outflux of x by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M13">View MathML</a>

(7)

The in silico labeling half-life of x is defined as the time when z drops to z0/2. We will show that this time equals the species half-life of x if its influx vin is zero. This property is independent from the amount of initially labeled entities, i.e. it holds for any z0/x0 ∈ ℝ+:

Proof:

Let x be determined by the processing with an unknown, potentially non-linear outflux vout and no influx vin = 0, i.e. v = vout,

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M14">View MathML</a>

Then, the kinetics of the label species z(t) is given by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M15">View MathML</a>

(8)

It can be shown that the factor <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M16">View MathML</a> is constant:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M17">View MathML</a>

Since this relation holds also true for t = 0, the proportionality constant is given by <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M18">View MathML</a> Then, equation (8) reads

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M19">View MathML</a>

Both processes x(t) and z(t) share the same half-life T1/2, since

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M20">View MathML</a>

This relation does not hold for processes with vin ≠ 0, because the fraction z/x becomes time-dependent as the labeling gets diluted, except for linear outfluxes as shown in the next section.

In silico labeling for linear processes

In this section, we prove that the label half-life coincides with the half-life of a species x which is produced by an unknown, potentially non-linear influx vin and is consumed by a linear process.

Proof:

Let <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M21">View MathML</a> be given by an unknown, potentially non-linear influx vin and a linear outflux, kx,

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M22">View MathML</a>

Then, the analytical half-life of x can be determined via

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M23">View MathML</a>

For the labeled system z it holds that

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M24">View MathML</a>

Creating the Extended Reaction Network

Some entities xi belong to the group of tracked, i.e. potentially labeled entities. Let us assume that they are given by x1, . . . , xα and untracked ones by xα+1, . . . , xm. Further, it can be assumed without loss of generality that (1) x1, . . . , xγ ≤ α are not complexes consisting of two or more tracked single entities, and (2) that the tracked single entities within each complex xγ+1, . . . , xα belong to the set x1, . . . , xγ. In the JAK-STAT example, S, pR_S, and pS belong to x1, . . . , xα and pS_pS to xα+1,. . . , xm as it contains two labeled single entities pS.

Creating additional entities xLF

A new set of labeled or free entities xLF is created based on the original x, by applying the following rules:

• Start with an empty set, xLF = {}

• Single entities: For each xi ∈{x1, . . . , xγ}, xLF is enlarged by a labeled <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M25">View MathML</a> and a free <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M26">View MathML</a> version of the original entity

• Complex entities: Each complex xi ∈ {xγ+1, . . . , xα} is decomposed into <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M27">View MathML</a> Due to the combinatorial multiplicity,

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M28">View MathML</a>

(9)

possible combinations using labeled <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M29','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M29">View MathML</a> and free <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M30">View MathML</a> entities are created, taking into account the order of the elements in the original complex xi, and are added to xLF. The complex pS_pS for instance leads to the four new complexes pSF_pSF, pSL_pSF, pSF_pSL, and pSL_pSL.

Creating additional reactions rLF

In order to create a new set of reactions rLF , the combinatorial multiplicity has to be applied not only to complexes but also to the ordered lists of reactants and products. Suppose an ordered list I of entities from the set {xi}1 ≤ iα with possible repetition, as for example the reactants of the reaction A + A + pA_pA A_A_pA_pA corresponding to I = (A, A, pA_pA). Summing up all single reactants and elements of the complexes leads to p single entities, in this case p = 4. Taking into account all combinations of labeled and free entities, 2p different lists can be derived, in the example

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M31">View MathML</a>

Without loss of generality, only the first δ reactions of the original system are assumed to affect a tracked entity. In these reactions, at least one reactant or product is a tracked entity. Then, a new set of reactions rLF can be established. Starting with the empty set rLF = {}, for each reaction ri {r1, . . . , rδ} with one or more reactants of tracked entities,

1. all reactants and products not belonging to the group of tracked entities are removed,

2. the combinatorial multiplicity approach is applied to the ordered list I of the remaining reactants leading to <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M32','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M32">View MathML</a>,

3. 2p reactions are added to rLF with reactants <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M32','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M32">View MathML</a> and the corresponding products.

4. the fluxes <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M33','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M33">View MathML</a> of the new reactions are given by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M34','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M34">View MathML</a>

(10)

Note that again the same symbol has been used for the entity name and its concentration. The sum over all weighting factors is 1.

Reactions ri ∈ {r1, . . . , rδ}without reactants produce only free entities, which simplifies the conversion of ri before adding to rLF: All untracked entities are removed, all xi are replaced by <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M26">View MathML</a>, and the flux is again given by equation (10).

When calculating the label half-life, products that coincide with the initially labeled entity are replaced by the corresponding free entity. This corresponds to removing the label and is necessary to avoid double-counting and to exclude upstream fluxes.

In order to calculate the label transit-time, entities entering the target pool must be released from their labeling, again, to avoid double-counting. Therefore, all labeled target entities are replaced in the reaction network rLF by their free counterparts. At the same time, a new product is added to those reactions where the target entity is a product to accumulate the returned label, RL.

Calculating the Label Half-Life and Transit-Time

Since the label half-life and transit-time characteristics are time-dependent, the label is not only injected at time-point 0, but the procedure is repeated for a series of time-points t (let xi be the source species):

1. Set all initial values for labeled entities and RL, if available, to 0. Set the initial value of free entities to the value of their counterpart in the original network.

2. Numerically integrate the ordinary differential equations corresponding to the extended reaction network {r, rLF} from 0 to t.

3. Apply a complete labeling of the source species: Set <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M35','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M35">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M36','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M36">View MathML</a> This step corresponds to the label injection.

4. Continue the numerical integration.

Threshold crossing at t" of the time-profiles <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M37','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M37">View MathML</a> and RL(t' > t) with <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M38','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M38">View MathML</a> defines the label half-life and label transit-time as t"-t, respectively. The threshold crossing is determined by linear interpolation of the discrete samples given by the numerical integration.

Profile Likelihood-based Confidence Intervals

We recently suggested a profile likelihood-based approach to determine simultaneous and separate confidence intervals for calibrated unknown model parameters [15]. In order to determine confidence intervals for the calculated label half-life and transit-times, the above procedure is not only repeated for a series of time-points, but also for a series of parameter settings. Each setting corresponds to one extreme point on the multi-dimensional manifold of acceptable parameter values, where one parameter has reached a lower or upper confidence threshold. By plotting all LHL or LTT profiles into one axis and creating an envelope between the largest and lowest values, a confidence interval for LHL and LTT is given.

Analytic half-lives for simple, isolated processes

The half-life T1/2(t) of simple and isolated biochemical reactions can be calculated analytically. Except for first-order processes, it usually depends on the concentration x0 = x(t0) at the time-point of interest t0 and is therefore time-dependent:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M39','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M39">View MathML</a>

(11)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M40','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M40">View MathML</a>

(12)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M41','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M41">View MathML</a>

(13)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M42','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M42">View MathML</a>

(14)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M43','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M43">View MathML</a>

(15)

The half-life calculation for a process of order n > 1 with <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M44">View MathML</a> is based on the integral form

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M45','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M45">View MathML</a>

(16)

In order to calculate the half-life for Michaelis-Menten kinetics, <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M46','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M46">View MathML</a> the following integral form is used which has been derived in [20], for x(t) with known x0 at t = t0:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M47','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M47">View MathML</a>

(17)

Panel A of Figure 1 displays the analytic results and their numerical approximation.

Results

In this section, the in silico labeling approach is applied to the JAK-STAT signaling pathway. The following mass action-based mechanistic model of the pathway has been calibrated to immunoblot measurements for Epo-stimulated BaF3-EpoR cells (model motivated by and data taken from [21]):

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M48','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M48">View MathML</a>

(18)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M49','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M49">View MathML</a>

(19)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M50','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M50">View MathML</a>

(20)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M51','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M51">View MathML</a>

(21)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M52','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M52">View MathML</a>

(22)

A smoothing spline approximation of the phosphorylated receptor served as the input function pR(t) triggering the phosphorylation of STAT (S pS). After dimerization (pS + pS pS-pS), the complexes enter the nucleus (pS_pS npS_npS). Then they dissociate and are dephosphorylated (npS_npS nS + nS). Finally, single STAT molecules leave the nucleus again (nS S). Model parameters were estimated using a Levenberg-Marquardt approach and the PottersWheel modeling software. The pools of total and phosphorylated cytoplasmic STAT have been used as observation functions. The kinetic parameters were estimated as k1 = 1.37, k2 = 0.22, k3 = 0.63, k4 = 0.59, and k5 = 0.59. The initial value of S was calibrated to 0.96 and the scaling factors for the observables to 1.45 for pS-obs and 0.98 for S-obs.

Labeled system

In order to determine the label half-life and transit-time of STAT, S is both, the initially labeled entity and the target pool. The flux of the label is illustrated in Figure 2. The time-courses of the original (solid blue) and labeled system (dashed red) are compared in Figure 3. In the beginning, both systems behave in the same manner. Then, the first wave of STAT molecules return from their cycle through the nucleus. Since they loose their label, the amount of labeled cytoplasmic STAT does not recover in contrast to the amount of STAT. After ~ 13 minutes, 50% of the initially labeled STAT molecules passed the nucleus at least once, as shown by the artificial pool of the returned label. The bimodal behavior of pSTAT exemplifies the first original signal wave and the secondary cycling effects. The in silico labeling approach allowed for discrimination between these two dynamics. In order to determine the transit-time for t > 0, the label is injected at a series of time-points, which is visualized in Figure 4.

thumbnailFigure 3. Time-courses of the original and labeled system. At t = 0, all STAT molecules are labeled and the original system (solid blue) coincides with the labeled one (dashed red). After the initial signal wave, STAT molecules return from the nucleus to the cytoplasm. Since the label transit-time for a complete cycle of cytoplasmic STAT is investigated, the returning molecules release their label, which is counted in an artificial entity (bottom right, green). Unlabeled molecules are depicted in yellow. Solid rose lines depict half labeled species, i.e. dimers of a labeled and an unlabeled molecule. They are shown only once, since the trajectory for example of pSL-pSF is identical to the one of pSF-pSL.

thumbnailFigure 4. Time-dependent label half-life and transit-time of cytoplasmic unphosphorylated STAT. A: In order to determine the label half-life and transit-time, a family of labeled STAT trajectories is calculated. For trajectory i, the amount of S(ti) label is injected into the system as SL at ti, with 0 ≤ ti ≤ 30. Therefore, the upper limit of subplot A matches the time-course of S for t ≤ 30 in Fig 3. B: Shown are the corresponding trajectories for returned label. The label half-life for trajectory i is the time period until <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M53','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M53">View MathML</a> drops below <a onClick="popup('http://www.biomedcentral.com/1752-0509/6/13/mathml/M54','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/6/13/mathml/M54">View MathML</a> C: The shortest label half-life is given as 0.6 ± 0.1 minutes. D: The label transit-time for a complete cycle of STAT molecules through the nucleus is the time-period until the returning label exceeds half of its initial value. Here, the minimum was estimated as 12 ± 2 minutes. Previously, it has been shown that the sojourn time of a single STAT5 molecule in the nucleus is about 6 minutes [21]. Our results indicate that on (median) average, a STAT5 molecule spends equal times in each compartment and requires about 12 min for one cycle from the cytoplasm to the nucleus and back. The 95% confidence intervals have been determined using the profile likelihood approach (PLE) and are displayed as grey envelope curves.

Label half-life and transit-time

Figure 4 depicts the label half-life of STAT as calculated from the time-course of SL. It reaches a minimum of 0.6 ± 0.1 minutes after ten minutes compared to a half-life of approximately 3 to 4 minutes at the starting point of the time-course analysis. For later time-points, the stimulus decreases (not shown) and STAT is no longer phosphorylated, resulting in an increased label half-life of STAT molecules. The minimum label transit-time for a complete cycle of STAT molecules was estimated as 12 ± 2 minutes.

Profile likelihood-based confidence intervals

In order to investigate how uncertainties in calibrated model parameter values propagate to the estimated time-characteristics, we applied the profile likelihood approach on an identifiable version of the model. The kinetic parameters involved in phosphorylation (k1), dimerization (k2), nuclear import (k3) and export (k5) were systematically varied consecutively within four orders of magnitude between 0.01kfit and 100kfit, with kfit being the parameter value for the best fit. For each variation, the other free parameters were calibrated resulting in a profile likelihood estimation (see Fig. S2 in additional file 1). All parameter settings corresponding to a crossing of the profile likelihood with the X2-threshold of the separate 95% confidence interval are used to recalculate the label half-life and transit-time. Figure 4C and 4D display the LHL and LTT 95% confidence interval by envelope curves. In case of the label half-life of cytoplasmic STAT, the confidence interval is very narrow allowing the LHL estimation within ± 0.1 minute for a range of label injection times between t = 0 and t = 20 minutes. The label transit-time has a wider confidence interval reflecting the larger number of reactions involved in a complete cycle of shuttling STAT.

Discussion and Conclusions

In this paper, the half-life of a species has been compared conceptually, analytically, and numerically to the half-life of a label in a hypothetical tracer experiment. Two time-characteristics, the label half-life and label transit time have been introduced, which capture the kinetics of a dynamical system on a higher level than e.g. single rate constants. Calculation of the time-characteristics and their profile likelihood-based confidence intervals for an identifiable pathway model showed that the approach is robust against parameter uncertainties. The quantities are calculated based on the novel in silico labeling method, which relies on an extended reaction network taking into account constraints concerning double-counting, upstream fluxes and combinatorial multiplicity. Our model-based in silico approach allows for insights into reaction networks that cannot be determined experimentally.

The proposed method provides important information for a wide spectrum of biological applications ranging from cell biology and pharmacokinetics to population dynamics. We applied it to a non-linear model of the cellular JAK-STAT signaling pathway, which allowed for calculating the time-dependent label half-life and transit-time of cytoplasmic STAT.

In summary our approach enables to calculate the amount of time a molecule spends in a certain state or compartment and therefore provides novel insights into the temporal scale of networks. This knowledge will have profound impact on drug design, as it offers the possibility to predict the life-time of a specific molecule and provides a basis to improve drug targeting.

Authors' contributions

TM and JB: Definition of the biological question, initiation, development, and implementation of the proposed method, writing of the manuscript. AR: Development of the proposed method, writing of the manuscript. SH, MS, SS, VB, UK, JT: Definition of the biological question, initiation of the method, critical discussion and contribution to manuscript. All authors read and approved the final manuscript.

Acknowledgements

This work was supported by the German Virtual Liver Network of the German Federal Ministry of Education and Research (BMBF, FKZ 0315766), the German Federal Ministry for Economy, the European Social Fund (BMWi, ESF, 03EGSBW-004), the Initiative and Networking Fund of the Helmholtz Association within the Helmholtz Alliance on Systems Biology/SBCancer Helmholtz, the NIH grant GM68762, and the German Federal Ministry for Education and Research (BMBF, FRISYS 0313921). This work was also supported by the Excellence Initiative of the German Federal and State Governments.

References

  1. Kitano H: Computational systems biology.

    Nature 2002, 420(6912):206-210.

    [PM:12432404]

    PubMed Abstract | Publisher Full Text OpenURL

  2. Hornberg JJ, Bruggeman FJ, Westerhoff HV, Lankelma J: Cancer: A Systems Biology disease. [http:/ / www.sciencedirect.com/ science/ article/ B6T2K-4J2TVW8-1/ 1/ a0e652495aabf78c62036fa1c08ab41a] webcite

    Biosystems 2006, 83:81-90. PubMed Abstract | Publisher Full Text OpenURL

  3. Paul Lee WN, Wahjudi PN, Xu J, Go VL: Tracer-based metabolomics: concepts and practices. [http://dx.doi.org/10.1016/j.clinbiochem.2010.07.027] webcite

    Clin Biochem 2010, 43(16-17):1269-1277. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  4. Bartelstone HJ, Mandel ID, Oshry E, Seidlin SM: Use of Radioactive Iodine as a Tracer in the Study of the Physiology of Teeth. [http://dx.doi.org/10.1126/science.106.2745.132-a] webcite

    Science 1947, 106(2745):132-133. PubMed Abstract | Publisher Full Text OpenURL

  5. Wang D, Shi J, Tan J, Jin X, Li Q, Kang H, Liu R, Jia B, Huang Y: Synthesis, characterization, and in vivo biodistribution of 125I-labeled Dex-g-PMAGGCONHTyr. [http://dx.doi.org/10.1021/bm200194s] webcite

    Biomacromolecules 2011, 12(5):1851-1859. PubMed Abstract | Publisher Full Text OpenURL

  6. Trembecka DO, Kuzak M, Dobrucki JW: Conditions for using FRAP as a quantitative technique-influence of the bleaching protocol. [http://dx.doi.org/10.1002/cyto.a.20866] webcite

    Cytometry A 2010, 77(4):366-370. PubMed Abstract | Publisher Full Text OpenURL

  7. Faeder JR, Blinov ML, Goldstein B, Hlavacek WS: Combinatorial complexity and dynamical restriction of network flows in signal transduction.

    Syst Biol (Stevenage) 2005, 2:5-15. Publisher Full Text OpenURL

  8. Golicnik M: Exact and approximate solutions for the decades-old Michaelis-Menten equation: Progress-curve analysis through integrated rate equations. [http://dx.doi.org/10.1002/bmb.20479] webcite

    Biochem Mol Biol Educ 2011, 39(2):117-125. PubMed Abstract | Publisher Full Text OpenURL

  9. Degen LP, Phillips SF: Variability of gastrointestinal transit in healthy women and men.

    Gut 1996, 39(2):299-305. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  10. Bergner PEE: [http://www.bergner.se/DMP/DMP%20Book%20edition%204.pdf] webcite

    A kinetics of macroscopic particles in open heterogeneous systems. 4th edition. Stockholm: Per-Erik E. Bergner; 2006. OpenURL

  11. Khinchin A: Mathematical Methods in the Theory of Queueing. Charles Griffin & Co. LTD, London; 1960. OpenURL

  12. Covell D, Berman M, Delisi C: Mean residence time - theoretical development, experimental determination, and practical use in tracer analysis.

    Mathematical Biosciences 1984, 72:213-244. Publisher Full Text OpenURL

  13. Weiss M: The relevance of residence time theory to pharmacokinetics.

    Eur J Clin Pharmacol 1992, 43(6):571-579. PubMed Abstract | Publisher Full Text OpenURL

  14. Veng-Pedersen P: Mean time parameters in pharmacokinetics. Definition, computation and clinical implications (Part II).

    Clin Pharmacokinet 1989, 17(6):424-440. PubMed Abstract | Publisher Full Text OpenURL

  15. Raue A, Kreutz C, Maiwald T, Bachmann J, Schilling M, Klingmüller U, Timmer J: Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood. [http://dx.doi.org/10.1093/bioinformatics/btp358] webcite

    Bioinformatics 2009, 25(15):1923-1929. PubMed Abstract | Publisher Full Text OpenURL

  16. Maiwald T, Timmer J: Dynamical Modeling and Multi-Experiment Fitting with PottersWheel.

    Bioinformatics 2008, 24(18):2037-2043. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  17. Becker V, Schilling M, Bachmann J, Baumann U, Raue A, Maiwald T, Timmer J, Klingmüller U: Covering a broad dynamic range: information processing at the erythropoietin receptor. [http://dx.doi.org/10.1126/science.1184913] webcite

    Science 2010, 328(5984):1404-1408. PubMed Abstract | Publisher Full Text OpenURL

  18. Kleiman LB, Maiwald T, Conzelmann H, Lauffenburger DA, Sorger PK: Rapid phospho-turnover by receptor tyrosine kinases impacts downstream signaling and drug binding. [http://dx.doi.org/10.1016/j.molcel.2011.07.014] webcite

    Mol Cell 2011, 43(5):723-737. PubMed Abstract | Publisher Full Text OpenURL

  19. Heinrich R, Schuster S: The Regulation of Cellular Systems. New York: Chapman & Hall; 1996. OpenURL

  20. Wagner JG: Properties of the Michaelis-Menten equation and its integrated form which are useful in pharmacokinetics.

    J Pharmacokinet Biopharm 1973, 1(2):103-121. PubMed Abstract | Publisher Full Text OpenURL

  21. Swameye I, Müller TG, Timmer J, Sandra O, Klingmüller U: Identification of nucleocytoplasmic cycling as a remote sensor in cellular signaling by databased modeling. [http://dx.doi.org/10.1073/pnas.0237333100] webcite

    Proc Natl Acad Sci USA 2003, 100(3):1028-1033. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL