Illustration of the method

Figure 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,

*S*^{L}/S^{F}, decreases in the cytoplasm. The total flux

of the first reaction, that is the phosphorylation of STAT molecules, can be divided into the flux

of labeled and the flux

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:

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

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

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

*r*_{j}, 1 ≤

*j *≤

*n*, corresponding to a set of coupled differential equations. The concentration change of each entity

*x*_{i}, 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]

Here,

*v*_{j }describes the flux of reaction

*j, a*_{ij }≥ 0 the stoichiometry of

*x*_{i }as a product in reaction

*j *and

*b*_{ij }≥ 0 the stoichiometry of

*x*_{i }as a reactant in reaction

*j*. We use the same symbol for an entity and its concentration, [

*x*_{i}]

*x*_{i}. The time-profile of each species can then be calculated for given initial values

and potentially driving input functions

*u*_{k}(

*t*). The flux

*v*_{j }of reaction

*j *may be a non-linear function of one or more species concentrations

*x*_{i }and externally defined

*u*_{k}. To improve readability, we omit explicitly denoting the time-dependency, i.e.

*x*_{i}(

*t*) is rather written as

*x*_{i}.

Analytical and numerical half-life calculation

The half-life of a species *x*_{i }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 *x*_{i},

with initial value

The whole system (4, 5) is solved either analytically or numerically and the species half-life of

*x*_{i }is given by

*T*_{1/2 }for which

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

and is therefore time-dependent. In this case, the above procedure is repeated for a series of different initial time-points

*t*_{0}. In a numerical integration, it is important to limit the maximum integrator step size for an accurate approximation of the

*y*^{0}/2 threshold crossing.

The half-life of a species *x*_{i }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 *x*_{i }has either no influx or when the outflux from *x*_{i }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 * {

*x*_{1}, . . . ,

*x*_{m}} is consumed only in one reaction. In

*in silico *labeling, the flux of the corresponding label

*z *depends on the outflux of

*x *by

The

*in silico *labeling half-life of

*x *is defined as the time when

*z *drops to

*z*_{0}*/*2. We will show that this time equals the species half-life of

*x *if its influx

*v*_{in }is zero. This property is independent from the amount of initially labeled entities, i.e. it holds for any

*z*_{0}*/x*_{0 } ^{+}:

Proof:

Let *x *be determined by the processing with an unknown, potentially non-linear outflux *v*_{out }and no influx *v*_{in }= 0, i.e. *v *= *v*_{out},

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

It can be shown that the factor

is constant:

Since this relation holds also true for

*t *= 0, the proportionality constant is given by

Then, equation (8) reads

Both processes *x*(*t*) and *z*(*t*) share the same half-life *T*_{1/2}, since

This relation does not hold for processes with *v*_{in }≠ 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 *v*_{in }and is consumed by a linear process.

Proof:

Let

be given by an unknown, potentially non-linear influx

*v*_{in }and a linear outflux,

*kx*,

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

For the labeled system *z *it holds that

Creating the Extended Reaction Network

Some entities *x*_{i }belong to the group of tracked, i.e. potentially labeled entities. Let us assume that they are given by *x*_{1}, . . . , *x*_{α }and untracked ones by *x*_{α+1}, . . . , *x*_{m}. Further, it can be assumed without loss of generality that (1) *x*_{1}, . . . , *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 *x*_{1}, . . . , *x*_{γ}. In the JAK-STAT example, *S, pR_S*, and *pS *belong to *x*_{1}, . . . , *x*_{α }and *pS*_*pS *to *x*_{α+1},. . . , *x*_{m }as it contains two labeled single entities *pS*.

Creating additional entities x^{LF}

A new set of labeled or free entities **x**^{LF }is created based on the original **x**, by applying the following rules:

• Start with an empty set, **x**^{LF }= {}

• Single entities: For each

*x*_{i }{

*x*_{1}, . . . ,

*x*_{γ}},

**x**^{LF }is enlarged by a labeled

and a free

version of the original entity

• Complex entities: Each complex

*x*_{i } {

*x*_{γ+1}, . . . ,

*x*_{α}} is decomposed into

Due to the combinatorial multiplicity,

possible combinations using labeled

and free

entities are created, taking into account the order of the elements in the original complex

*x*_{i}, and are added to

**x**^{LF}. The complex

*pS*_{_}*pS *for instance leads to the four new complexes

*pS*^{F}_pS^{F}, pS^{L}_pS^{F}, pS^{F}_pS^{L}, and

*pS*^{L}_pS^{L}.

Creating additional reactions r^{LF}

In order to create a new set of reactions **r**^{LF }, 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 {*x*_{i}}_{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, 2^{p }different lists can be derived, in the example

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

**r**^{LF }can be established. Starting with the empty set

**r**^{LF }= {}, for each reaction

*r*_{i } **{***r*_{1}, . . . ,

*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

,

3. 2

^{p }reactions are added to

**r**^{LF }with reactants

and the corresponding products.

4. the fluxes

of the new reactions are given by

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

*r*_{i } {

*r*_{1}, . . . ,

*r*_{δ}}without reactants produce only free entities, which simplifies the conversion of

*r*_{i }before adding to

**r**^{LF}: All untracked entities are removed, all

*x*_{i }are replaced by

, 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 **r**^{LF }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 *x*_{i }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**, **r**^{LF}} from 0 to *t*.

3. Apply a complete labeling of the source species: Set

and

This step corresponds to the label

*injection*.

4. Continue the numerical integration.

Threshold crossing at

*t*" of the time-profiles

and RL(

*t' **> t*) with

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 *T*_{1/2}(*t*) of simple and isolated biochemical reactions can be calculated analytically. Except for first-order processes, it usually depends on the concentration *x*_{0 }= *x*(*t*_{0}) at the time-point of interest *t*_{0 }and is therefore time-dependent:

The half-life calculation for a process of order

*n >*1 with

is based on the integral form

In order to calculate the half-life for Michaelis-Menten kinetics,

the following integral form is used which has been derived in [

20], for

*x*(

*t*) with known

*x*_{0 }at

*t *=

*t*_{0}:

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