In silico labeling reveals the timedependent label halflife and transittime in dynamical systems
 Thomas Maiwald†^{1, 2}Email author,
 Julie Blumberg†^{3, 4},
 Andreas Raue^{1},
 Stefan Hengl^{1},
 Marcel Schilling^{5},
 Sherwin KB Sy^{6},
 Verena Becker^{2, 5, 7},
 Ursula Klingmüller^{5, 7} and
 Jens Timmer^{1, 8, 9, 10}
DOI: 10.1186/17520509613
© Maiwald et al.; licensee BioMed Central Ltd. 2012
Received: 10 May 2011
Accepted: 27 February 2012
Published: 27 February 2012
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 timeperiod a protein is accessible to other molecules in a certain state  its halflife  and (2) the time it spends when passing through a subsystem  its transittime. We discuss two approaches to quantify the halflife, present the novel method of in silico labeling, and introduce the label halflife and label transittime. 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 doublecounting 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 halflife and transittime.
Results
Application to the JAKSTAT signaling pathway in Epostimulated BaF3EpoR cells enabled the calculation of the timedependent label halflife and transittime of STAT species. The results were robust against parameter uncertainties.
Conclusions
Our approach renders possible the estimation of species and label halflives and transittimes. It is applicable to large nonlinear systems and an implementation is provided within the PottersWheel modeling framework (http://www.potterswheel.de).
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 halflife, applies the halflife concept not to a species, but to a virtual label attached to the species. The second one, the label transittime, is the timeperiod 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 HalfLife
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 iodine125 [4, 5] or green fluorescent proteintagged 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 deadends, avoid doublecounting of cycling label, and to restrict the label to a subnetwork 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].
While for simple systems the species halflife can be determined analytically, the symbolic integration of a MichaelisMenten 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 halflife.
Label TransitTime
Transittimes 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 transittime denotes the time spent on average in a subsystem [10], while the mean sojourntime also takes into account the probability that the subsystem is entered at all [11]. In pharmacokinetics, the socalled 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 concentrationtime profile of a drug [14].
We here introduce the label transittime (LTT) from a source to a target pool in a chemical reaction network as the timeperiod 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 halflife, it is important to distinguish entities residing in the source pool at t = 0 from other entities entering the source pool at later timepoints. When calculating transittimes, 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 halflife 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 doublecounting the labeled species. In the case of transittime 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 bookkeeping 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 Likelihoodbased Confidence Intervals
Recently, we suggested a profile likelihoodbased 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 timedependent label halflife and transittime characteristics.
Implementation
All concepts have been implemented within the PottersWheel modeling and parameter estimation framework that is available from http://www.potterswheel.de[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.
In the next section, the proposed labeling method is illustrated for the JAKSTAT signal transduction pathway and afterwards described in detail. After proving the equality of species and label halflife for isolated or linear processes, a fitted model of the JAKSTAT pathway is used to determine the label halflife of unphosphorylated STAT and its label transittime when cycling through the nucleus of a cell.
Methods
Illustration of the method
This procedure is repeated for a series of timepoints t in order to determine LHL(t) and LTT(t) for all timepoint of interest.
Terminology
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 timeprofile of each species can then be calculated for given initial values ${x}_{i}^{0}={x}_{i}\left({t}_{0}\right)$ and potentially driving input functions u_{ k } (t). The flux v_{ j } of reaction j may be a nonlinear function of one or more species concentrations x_{ i } and externally defined u_{ k } . To improve readability, we omit explicitly denoting the timedependency, i.e. x_{ i } (t) is rather written as x_{ i } .
Analytical and numerical halflife calculation
Note that a halflife 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 halflife. Otherwise, the halflife depends on the initial concentration ${x}_{i}^{0}$ and is therefore timedependent. In this case, the above procedure is repeated for a series of different initial timepoints 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 halflife 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 halflife as a means to determine a timecharacteristic which is motivated by laboratory tracer experiment with the additional property to avoid tracerdouble counting in kinetic cycles.
In silico labeling halflife for isolated processes
The in silico labeling halflife of x is defined as the time when z drops to z_{0}/2. We will show that this time equals the species halflife 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:
This relation does not hold for processes with v_{ in } ≠ 0, because the fraction z/x becomes timedependent 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 halflife coincides with the halflife of a species x which is produced by an unknown, potentially nonlinear influx v_{ in } and is consumed by a linear process.
Proof:
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 JAKSTAT 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 ${x}_{i}^{L}$ and a free ${x}_{i}^{F}$ version of the original entity

Complex entities: Each complex x_{ i }∈ {x_{γ+1}, . . . , x_{α}} is decomposed into ${n}_{1}^{i}{x}_{1},\dots ,{n}_{\gamma}^{i}{x}_{\gamma}.$ Due to the combinatorial multiplicity,${2}^{{\displaystyle {\sum}_{j=1}^{\gamma}{n}_{j}^{i}}}$(9)
possible combinations using labeled ${x}_{j}^{L}$ and free ${x}_{j}^{F}$ 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 }
 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 ${I}_{1},...,{I}_{{2}^{p}}$,
 3.
2 ^{ p } reactions are added to r ^{ LF } with reactants ${I}_{1},...,{I}_{{2}^{p}}$ and the corresponding products.
 4.the fluxes ${v}_{1}^{LF},\dots ,{v}_{j}^{LF}$ of the new reactions are given by${v}_{j}^{LF}=\frac{{\prod}_{{x}_{k}^{LF}\in {I}_{j}}{x}_{k}^{LF}}{{\prod}_{{x}_{k}\in I}{x}_{k}}\phantom{\rule{2.22198pt}{0ex}}{v}_{i},\phantom{\rule{2.22198pt}{0ex}}\phantom{\rule{2.22198pt}{0ex}}1\le j\le {2}^{p}.$(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 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 ${x}_{i}^{F}$, and the flux is again given by equation (10).
When calculating the label halflife, 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 doublecounting and to exclude upstream fluxes.
In order to calculate the label transittime, entities entering the target pool must be released from their labeling, again, to avoid doublecounting. 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 HalfLife and TransitTime
 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 ${x}_{i}^{L}\left(t\right)={x}_{i}\left(t\right)$ and ${x}_{i}^{F}\left(t\right)=0.$ This step corresponds to the label injection.
 4.
Continue the numerical integration.
Threshold crossing at t" of the timeprofiles ${x}_{i}^{L}\left({t}^{\prime}>t\right)$ and RL(t' > t) with ${x}_{i}^{L}\left(t\right)/2$ defines the label halflife and label transittime as t"t, respectively. The threshold crossing is determined by linear interpolation of the discrete samples given by the numerical integration.
Profile Likelihoodbased Confidence Intervals
We recently suggested a profile likelihoodbased approach to determine simultaneous and separate confidence intervals for calibrated unknown model parameters [15]. In order to determine confidence intervals for the calculated label halflife and transittimes, the above procedure is not only repeated for a series of timepoints, but also for a series of parameter settings. Each setting corresponds to one extreme point on the multidimensional 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 halflives for simple, isolated processes
Panel A of Figure 1 displays the analytic results and their numerical approximation.
Results
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 LevenbergMarquardt 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 k_{1} = 1.37, k_{2} = 0.22, k_{>3} = 0.63, k_{4} = 0.59, and k_{5} = 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
Label halflife and transittime
Figure 4 depicts the label halflife of STAT as calculated from the timecourse of S^{ L } . It reaches a minimum of 0.6 ± 0.1 minutes after ten minutes compared to a halflife of approximately 3 to 4 minutes at the starting point of the timecourse analysis. For later timepoints, the stimulus decreases (not shown) and STAT is no longer phosphorylated, resulting in an increased label halflife of STAT molecules. The minimum label transittime for a complete cycle of STAT molecules was estimated as 12 ± 2 minutes.
Profile likelihoodbased confidence intervals
In order to investigate how uncertainties in calibrated model parameter values propagate to the estimated timecharacteristics, we applied the profile likelihood approach on an identifiable version of the model. The kinetic parameters involved in phosphorylation (k_{1}), dimerization (k_{2}), nuclear import (k_{3}) and export (k_{5}) were systematically varied consecutively within four orders of magnitude between 0.01k_{ fit } and 100k_{ fit } , with k_{ fit } 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 X^{2}threshold of the separate 95% confidence interval are used to recalculate the label halflife and transittime. Figure 4C and 4D display the LHL and LTT 95% confidence interval by envelope curves. In case of the label halflife 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 transittime 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 halflife of a species has been compared conceptually, analytically, and numerically to the halflife of a label in a hypothetical tracer experiment. Two timecharacteristics, the label halflife 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 timecharacteristics and their profile likelihoodbased 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 doublecounting, upstream fluxes and combinatorial multiplicity. Our modelbased 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 nonlinear model of the cellular JAKSTAT signaling pathway, which allowed for calculating the timedependent label halflife and transittime 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 lifetime of a specific molecule and provides a basis to improve drug targeting.
Notes
Declarations
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, 03EGSBW004), 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.
Authors’ Affiliations
References
 Kitano H: Computational systems biology. Nature. 2002, 420 (6912): 206210. 10.1038/nature01254. [PM:12432404]View Article
 Hornberg JJ, Bruggeman FJ, Westerhoff HV, Lankelma J: Cancer: A Systems Biology disease. Biosystems. 2006, 83: 8190. 10.1016/j.biosystems.2005.05.014. [http://www.sciencedirect.com/science/article/B6T2K4J2TVW81/1/a0e652495aabf78c62036fa1c08ab41a] 10.1016/j.biosystems.2005.05.014View Article
 Paul Lee WN, Wahjudi PN, Xu J, Go VL: Tracerbased metabolomics: concepts and practices. Clin Biochem. 2010, 43 (1617): 12691277. 10.1016/j.clinbiochem.2010.07.027. [http://dx.doi.org/10.1016/j.clinbiochem.2010.07.027] 10.1016/j.clinbiochem.2010.07.027View Article
 Bartelstone HJ, Mandel ID, Oshry E, Seidlin SM: Use of Radioactive Iodine as a Tracer in the Study of the Physiology of Teeth. Science. 1947, 106 (2745): 132133. [http://dx.doi.org/10.1126/science.106.2745.132a]View Article
 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 125Ilabeled DexgPMAGGCONHTyr. Biomacromolecules. 2011, 12 (5): 18511859. 10.1021/bm200194s. [http://dx.doi.org/10.1021/bm200194s] 10.1021/bm200194sView Article
 Trembecka DO, Kuzak M, Dobrucki JW: Conditions for using FRAP as a quantitative techniqueinfluence of the bleaching protocol. Cytometry A. 2010, 77 (4): 366370. [http://dx.doi.org/10.1002/cyto.a.20866]View Article
 Faeder JR, Blinov ML, Goldstein B, Hlavacek WS: Combinatorial complexity and dynamical restriction of network flows in signal transduction. Syst Biol (Stevenage). 2005, 2: 515. 10.1049/sb:20045031.View Article
 Golicnik M: Exact and approximate solutions for the decadesold MichaelisMenten equation: Progresscurve analysis through integrated rate equations. Biochem Mol Biol Educ. 2011, 39 (2): 117125. 10.1002/bmb.20479. [http://dx.doi.org/10.1002/bmb.20479] 10.1002/bmb.20479View Article
 Degen LP, Phillips SF: Variability of gastrointestinal transit in healthy women and men. Gut. 1996, 39 (2): 299305. 10.1136/gut.39.2.299.View Article
 Bergner PEE: A kinetics of macroscopic particles in open heterogeneous systems. 2006, Stockholm: PerErik E. Bergner, [http://www.bergner.se/DMP/DMP%20Book%20edition%204.pdf]4
 Khinchin A: Mathematical Methods in the Theory of Queueing. 1960, Charles Griffin & Co. LTD, London
 Covell D, Berman M, Delisi C: Mean residence time  theoretical development, experimental determination, and practical use in tracer analysis. Mathematical Biosciences. 1984, 72: 213244. 10.1016/00255564(84)901111.View Article
 Weiss M: The relevance of residence time theory to pharmacokinetics. Eur J Clin Pharmacol. 1992, 43 (6): 571579. 10.1007/BF02284953.View Article
 VengPedersen P: Mean time parameters in pharmacokinetics. Definition, computation and clinical implications (Part II). Clin Pharmacokinet. 1989, 17 (6): 424440. 10.2165/0000308819891706000005.View Article
 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. Bioinformatics. 2009, 25 (15): 19231929. 10.1093/bioinformatics/btp358. [http://dx.doi.org/10.1093/bioinformatics/btp358] 10.1093/bioinformatics/btp358View Article
 Maiwald T, Timmer J: Dynamical Modeling and MultiExperiment Fitting with PottersWheel. Bioinformatics. 2008, 24 (18): 20372043. 10.1093/bioinformatics/btn350.View Article
 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. Science. 2010, 328 (5984): 14041408. 10.1126/science.1184913. [http://dx.doi.org/10.1126/science.1184913] 10.1126/science.1184913View Article
 Kleiman LB, Maiwald T, Conzelmann H, Lauffenburger DA, Sorger PK: Rapid phosphoturnover by receptor tyrosine kinases impacts downstream signaling and drug binding. Mol Cell. 2011, 43 (5): 723737. 10.1016/j.molcel.2011.07.014. [http://dx.doi.org/10.1016/j.molcel.2011.07.014] 10.1016/j.molcel.2011.07.014View Article
 Heinrich R, Schuster S: The Regulation of Cellular Systems. 1996, New York: Chapman & HallView Article
 Wagner JG: Properties of the MichaelisMenten equation and its integrated form which are useful in pharmacokinetics. J Pharmacokinet Biopharm. 1973, 1 (2): 103121. 10.1007/BF01059625.View Article
 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. Proc Natl Acad Sci USA. 2003, 100 (3): 10281033. 10.1073/pnas.0237333100. [http://dx.doi.org/10.1073/pnas.0237333100] 10.1073/pnas.0237333100View Article
Copyright
This article is published under license to 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.