Qi X. A Stochastic Version of Corticosteriod Pharmacogenomic Model. AAPS Journal. 2005; 7(1): E134-E140. DOI:  10.1208/aapsj070113

Correspondence to:
Xiaohong Qi
Tel: (954) 358-6160
Fax: (954) 358-6352
Email: xiaohong.qi@andrx.com

Received: May 28, 2004;  Accepted: November 15, 2004;  Published: June 23, 2005

Abstract

The purpose of this study was to develop a stochastic version of corticosteriod fifth generation pharmacogenomic model. The Gillespie algorithm was used to generate the independent time courses of the receptor messenger RNA (mRNA). Initial parameters for the stochastic simulation were adapted from the study by Jin et al. The result obtained from the proposed stochastic model showed an overall agreement with the deterministic fifth generation model. This study suggested that because the stochastic model takes into account the “noise” nature of gene regulation, it would have potential application in pharmacogenomic modeling.

Keywords: pharmacogenomic modeling, corticosteroid, stochastic, gene, Gillespie algorithm, Monte Carlo simulation

Introduction

Several mathematical models for corticosteriod (CS) pharmacokinetics-pharmacodynamics-pharmacogenomics have been proposed on the basis of the experimental observations.^1-7 These models are of a deterministic nature and take the form of a system of the coupled ordinary differential equations. The CS gene regulatory networks and cell signaling pathways exhibited nonlinear behavior, reflected by the complex mechanisms of interactions, such as enzymatic actions, gene expression, and various feedback loops. The gene network with a high nonlinearity could possess multiple stable states and bifurcations.⁸ Meanwhile, the numbers of molecules involved in the CS gene expression are low, for example, the baseline messenger RNA (mRNA) is 25.8 fmol/g liver, and the baseline glucocorticoid receptor is 540.7 fmol/mg protein. Based on the recent experimental and theoretical studies, the stochastic nature of the biochemical reactions must lead to large fluctuations.^9-10 Studies also confirmed that the phenotypic variability of organisms is highly related to the inherent stochasticity that operates at a basal level of gene expression.^11-14 The level of expression from the same gene could vary significantly from one cell to another within a genetically identical colony.^15-16 Furthermore, numerous findings showed that the origin of such variability among isogenic population is largely attributed to stochastic phenomena. The stochastic nature of the gene network motivated the development of a stochastic version of CS pharmacogenomics.

While the importance of fluctuations in the gene expression was stressed over 30 years ago,¹⁷ the recent expanding experimental evidence has revived interest in the use of stochastic simulation techniques to model gene regulatory networks.^14,18-20 Simulations algorithms are critical to successfully model the biological systems. Swain et al¹⁹ proposed a theoretical framework that enables interoperation of experimental measurements of stochasticity in gene expression. Their approach may have potential application in pharmacogenomic study. Stochastic fluctuations in pharmacokinetics and pharmacodynamics have drawn the attention of several research groups.^21-24 The use of the stochastic approach to address “macro” and “observable” variability problems in population pharmacokinetics and pharmacodynamics has been well demonstrated.^22-23 However, few reports use the stochastic approach to investigate the “micro” variability in pharmacogenomics, probably owing to a limited data resource.

When fluctuations arise from the small number of reactant molecules, the stochastic simulation algorithm developed by Gillespie is considered the gold standard for modeling.^25-26 The advantage of this algorithm is that it generates an ensemble of trajectories with correct statistics for a set of biochemical reactions.

The purpose of this study is to develop a stochastic version of the CS fifth generation pharmacogenomic model. The study is broadly divided into 3 sections. The first section introduces the Gillespie algorithm. In addition, a comparison of the stochastic approach with the deterministic approach was made in order to appreciate the relative merit of the stochastic method in the analysis of gene regulation data. The second section presents the stochastic results of CS pharmacogenomics. The third section discusses the effect of stochasticity on pharmacogenomics and the potential application of the stochastic approach in pharmacogenomics research.

Background of the Gillespie Algorithm

The Gillespie algorithm has been used to determine the dynamics of the biochemical or genetic systems in the presence of molecular noises.^9,11 Since the CS fifth generation pharmacogenomic model is so complex, a simple example was first used to demonstrate the feature of Gillespie algorithm.

Consider the following reactions, in which E represents the enzyme, S represents the substrate, ES represents the intermediate complex, and P represents the final product:

Reaction 1: E + S→ ES
Reaction 2: ES→ E + S
Reaction 3: ES→ E + P

The initial numbers of molecules for E, S, ES, and P are [E]₀, [S]₀, [ES]₀, and [P]₀, respectively. The simulation time was set to T. The reaction rate constants for Reactions 1 to 3 are k₁, k₂, and k₃, respectively.

In the Gillespie algorithm, the state of the system is updated by determining (1) the time to the next reaction, and (2) which reaction will occur next. A single iteration of Gillespie algorithm for the above reactions was performed as follows:

Step 1

Generate 2 random numbers U₁ and U₂ uniformly distributed over (0,1). U₁ determines which reaction in the system is going to occur; U₂ determines the waiting time for this reaction.

Step 2

Execute one elementary Reaction 2. Set the number of ES to be [ES]₀ − 1; set the number of P to be [P]₀ + 1; and set the number of E to be [E]₀ + 1.

Set simulation time to T + t.

Step 3

The probability of choosing a given reaction is proportional to the product of its rate constant and numbers of substrate molecules. The waiting time for the reaction is calculated as a random variable generated according to exponential distribution with parameter a₀: t = ( 1 / a 0 ) ln ( 1 / U 2 ) . U₂ is uniformly distributed over (0,1), and a₀ = k₁× [E] × [S] + k₂ × [ES] + k₃ × [ES].

We now apply Gillespie algorithm to N chemical species that react through J reactions. X_j is the number of molecules of type j in the system. For each reaction we assumed that there existed a rate constant k_μ. The propensity of each reaction is expressed as a μ = k μ × H μ , where H_μ represents the possible combination of reactants. The first question is addressed by considering the distribution of the time intervals τ until the next reaction, P(τ). The mathematical equation for P(τ) is as follows:

P ( τ ) = a μ exp [ - τ ∑ μ a μ ]

(1)

The choice of the next reaction is made based upon the discrete distribution as follows:

P ( μ = μ ' ) = a μ ' / ∑ μ a μ

(2)

Table 1 summarized the procedures to perform Gillespie algorithm.

In order to appreciate the merit of the stochastic method in the analysis of gene regulation data, we compared it with the deterministic approach. The deterministic approach is based on the assumption that the concentration varies continuously and continually over time. In other words, the concentrations are real numbers. The variation of concentrations shows no discontinuities or “jumps.” This assumption is valid as long as the concentrations are high, while it is invalid in the cellular-based systems. For example, a single DNA molecule, RNA, and protein molecules ranged from tens to a few hundred. Besides, studies showed that the mRNA production is quantal and is produced in random pulses.²⁷ As a result, the discrete nature of change needs to be captured at a low level condition. In contrast to the deterministic approach, the stochastic simulations consider that the number of reactants is discrete and protein production occurs in short “bursts” at random time intervals rather than in a continuous manner.²⁸

In a deterministic setting, an initial condition would always guide a cell through a particular trajectory with no scope for flexibility of response. However, experiments showed that the gene expression variability and the accompanying noise make the cell flexible and help it to adapt to varying environmental conditions. Elowitz and coworkers²⁹ constructed strains of Escherichia coli for detecting noise and discriminating between the intrinsic and extrinsic contributors of noise. They found that both stochasticity, inherent in the biochemical process of gene expression (intrinsic noise), and fluctuations in other cellular components (extrinsic noise) contribute substantially to overall variation. Factors such as transcription rate, regulatory dynamics, and genetic factors control the amplitude of noise.

Identical initial conditions such as concentrations of chemical species, temperature, pressure, and so forth, can also produce qualitatively different outcomes in the temporal evolution of a regulatory network. When 2 independently produced regulatory proteins acting at low cellular concentrations competitively control a switch point in a pathway, stochastic variations in their concentrations can produce probabilistic pathway selection, so that an initially homogeneous cell population partitions into distinct phenotypic subpopulations.³⁰ The role of the noise must be taken into consideration in gene regulation. Levin³¹ showed that noise is a source of nongenetic individuality in the chemotactic response of E coli. This study is a good example of the importance and potential application of the stochastic approach in pharmacogenomics, as the interaction of drug with the individual cell will influence the efficacy of drug.

Overall, stochastic modeling takes a mesoscopic view of the system and keeps track of every molecule in the system and is close to real world accuracy. In the next section, the Gillespie algorithm is used to gather a stochastic version of CS pharmacogenomics.

Implementation

The Fifth Generation Pharmacogenomic Model of Corticosteroid

Figure 1 shows the known fifth generation model.^5-7 Here it is used without any modification. The fifth generation model represented a prototype for the molecular mechanism of CS that was based on the negative feedback exerted by a drug-receptor on the expression of its gene. Briefly, CS and cytosolic glucocorticoid receptor (GR) formed a CS-receptor complex (DR). DR was translocated into the nucleus and formed a drug-receptor in the nucleus (DR[N]). DR(N) feed back to inhibit the production of mRNA. k_on is the second-order rate constant of the formation of DR. k_T is the translocation rate. IC_{50_Rm} is the concentration of DR(N) at which the synthesis rate of mRNA drops to 50% of its baseline value. k_{s_Rm} is the zero-order rate of mRNA synthesis. k_{d_Rm} is the first-order rate of mRNA. k_{s_R} is the receptor synthesis rate constant. k_{d_R} is the receptor degradation rate constant. k_re is the DR(N) degradation rate constant. R_f is the fraction of free receptor being recycled.

Figure 1.. Fifth generation CS model.5-7 The dotted line represented inhibition of gene transcription by the DR(N) via an indirect mechanism.

The Stochastic Version of Corticosteroid Pharmacogenomics

Table 2 shows the stochastic version of CS pharmacogenomics. The first column is the reaction number. The second column is the reaction sequences. The propensity of each reaction is given in the third column. The propensity is a product of the reaction rate and a function of the number of the molecules involved. The last column indicates the changes in the numbers of molecules taking part in the different reactions. The sequence of reactions corresponds to the mechanism underlying the fifth generation model. In the 2-compartment model, k₁₂, k₂₁, and k₁₀ are the distribution rate constants.^5-7

A total of 11 types of biochemical reactions occurred in the CS gene network, with 5 species involved. According to the Gillespie algorithm, at each time step the algorithm randomly determines the reaction that takes place according to its probability, as well as the time interval to the next reaction step. The Gillespie algorithm was implemented in the MATLAB (MATLAB Version 6.5, The Mathworks Inc, Natick, MA). Monte Carlo simulations were run on a Compaq 866 (Hewlett Packard Ltd) and the program took less than 15 minutes to carry out. Initial parameters for the simulations were adapted from Jin et al,⁷ in which a similar system to Ramakrishnan et al^5,6 was modeled: k_on = 0.00329 nM ⁻¹ h ⁻¹; k_T = 0.63 h ⁻¹; IC_{50_Rm} = 26.2 fmol/mg protein; k_{s_Rm} = 2.9 fmol/g liver/h; k_{d_R} = 0.0572 h ⁻¹; k_re = 0.57 h ⁻¹; R_f = 0.49. The constants k_{d_Rm} and k_{s_R} were calculated according to the baseline equations, Equations 5 to 7 (equations not shown). Our direct reimplementation of the fifth generation model yielded similar results (data not shown).

Results

As shown in Table 2, we attributed to each linear or nonlinear term of the kinetic equations a probability of the occurrence of the corresponding reaction. The numbers of molecules of the different reacting species as well as the probabilities are updated at each time step. For example, Reaction 1 corresponds to the transcription of gene into mRNA. If this reaction occurs, the number of mRNA will increase by 1. Reaction 3 and Reaction 8 both contributed to the increase of the number of receptor protein by 1, but each has different probability.

Reaction 5 results in increasing by 1 the number of DR and decreasing by 1 the number of receptor (denoted as R). D_T in Table 2 represents the drug level in the peripheral compartment. In the 2-compartment model, k₁₂, k₂₁, and k₁₀ are the distribution rate constants.^5-7

One of the features of the Gillespie algorithm is used to generate many independent time courses for the system of coupled chemical reactions under investigation.

Independent trajectories may be used to study distribution of various system parameters.^25-26 Figure 2 illustrates the simulated mRNA course using the stochastic model. The data presented in Figure 2 show that this down-regulation pattern was recovered when using the stochastic model (ie, upon intravenous administration of corticosteriod drug, the receptor mRNA rapidly decline). The minimum mRNA value ranged from 5 to ∽25 nM with a median range from 10 to ∽15 nM. After 24 hours, ∽80% of mRNA returned to 20 to ∽30 nM.

Figure 2..  (A) The simulated temporal receptor mRNA course from the fifth generation CS model; (B) A histogram of minimum molecular number of mRNA; (C) A histogram of mRNA after 24 hours. Note: The histograms in panel B and panel C were constructed from frequency. The intervals of molecular number of mRNA are shown on the x-axis and the number of scores in each interval is represented by the height of a rectangle located above the interval.

Discussion

The inherent stochasticity of biochemical processes such as transcription and translation generates “intrinsic” noise; in addition, fluctuations in the amounts or states of other cellular components lead indirectly to variations in the expression of a particular gene and thus represent “extrinsic” noise.^9-11,13,15 The stochastic approach would be a favorable tool in the understanding of cell kinetics/dynamics and the analysis of the drug resistance during long-term drug therapy. For example, HIV is known to be a highly variable virus, with errors in reverse transcription during cell infection frequently resulting in the progeny virus differing from the parent.³² Specific viral populations sometimes have different susceptibility to antiretroviral drugs.³³ Stochastic processes showed a strong influence on the HIV-1 evolution during suboptimal-inhibitor therapy.³⁴ Random chance could play a crucial role in determining the existence of viral subpopulations having more substitutions and in people on therapy, in which there is a low new cell infection rate. Phillips et al³⁵ developed a stochastic model that mimics the dynamic processes thought to be operating in an infected individual who has been treated with antiretroviral drugs. In their model, a random number generator was used to determine whether mutations occur in any one round of replication and to sample the population created in the next generation. Their model illustrated that chance alone could result in a completely different outcome for patients who started off in identical situations.

An individual gene has the capability to adapt to varying environments. It has been postulated that through the mechanism of regulatory feedback loops, networks can impose certain average behavior even in the presence of randomness. The simulation of a molecular model of the circadian rhythms of the PER protein in Drosophila and of the FRQ protein in Neurospora illustrates that robust circadian oscillations can occur even with a limited number of mRNA and protein molecules (ie, in the range of tens and hundreds, respectively).³⁶ This finding has an implication in pharmacogenomics studies. Since most pharmacological effects reflect interactions of drug with populations of cell, when the concentrations of reactants are high, the fluctuations in the cell can be averaged out over the whole population, leaving the deterministic formulations intact. However, caution must be taken after intravenous administration of short half-life drug, in which case drug concentration will decay exponentially, and the intrinsic and extrinsic variability of the gene will take place. In a recent study, Gillespie³⁷ presented ideas to link stochastic and deterministic simulations, which might have potential application to deal with this problem.

Investigations of the fluctuations in single cells would play an important role in finding suitable drug candidates to solve biological rhythm problems (eg, sleep disorder, hormone secretion dysfunction). Sleep, hormonal secretions, and body temperature are body functions in a 24-hour cycle. Recent molecular dissections of the circadian biological clock system have revealed that oscillation in the transcription of specific clock genes plays a central role in the generation of circadian rhythms.^38-40 Several drugs, for example, interferon (IFN), can affect the expression of clock genes, resulting in alteration of the 24-hour rhythms in both physiology and behavior. Koyanagi et al³⁸ found that expression of vascular endothelial growth factor (VEGF) in hypoxic tumor cells was affected by the circadian organization of molecular clockwork. Because angiogenesis is essential for tumor growth and metastasis, inhibition of angiogenesis has emerged as a new therapy to treat cancers. The core circadian oscillator is composed of an autoregulatory transcription-translation feedback loop in which CLOCK and BMAL1 are positive regulators, and Period and Cryptochrome genes act as negative ones. The levels of VEGF mRNA in tumor cells implanted in mice rose substantially in response to hypoxia, but the levels fluctuated rhythmically in a circadian fashion. Luciferase reporter gene analysis revealed that Period2 and Cryptochrome1, whose expression in the implanted tumor cells showed a circadian oscillation, inhibited the hypoxia-induced VEGF promoter activity. Their study suggests that the negative limbs of the molecular loop periodically inhibit the hypoxic induction of VEGF transcription, resulting in the circadian fluctuation of the mRNA expression.

The stochastic models could provide probabilistic distributions of gene expression levels. This property can be used for describing the variation of expression products from cell to cell, which gives another possible means to investigate the so-called “random error” in gene expression data, especially for probes with low intensity levels.⁴¹ The reliability of data mining results can be assessed by comparing random errors estimated from observed data with the variation predicted by stochastic models.⁴²

Conclusion

The application of the stochastic approach in CS pharmacogenomic modeling was demonstrated in this study. Monte Carlo simulation using the Gillespie algorithm resulted in an overall agreement with the deterministic fifth generation model. This study suggests that because the stochastic model takes into account the “noise” nature of gene regulation, it would have potential application in pharmacogenomic modeling.

Acknowledgments

The author thanks Dr Huimin Hou (NPERC, China) for encouragement and Dr Jin Yan Jin (SUNY, Buffalo, NY) for reading the first draft of this manuscript and sending me a copy of her doctoral thesis. The author is grateful for Prof Ho-Leung Fung (SUNY) and for journal reviewers whose comments and suggestions substantially improved the quality of this manuscript.

References

1. Jusko WJ.  Receptor-mediated pharmacodynamics of corticosteroids. Prog Clin Biol Res. 1994;387:261-270.
PubMed 

2. Sun YN, DuBois DC, Almon RR, Jusko WJ.  Fourth-generation model for corticosteroid pharmacodynamics: a model for methylprednisolone effects on receptor/gene-mediated glucocorticoid receptor down-regulation and tyrosine aminotransferase induction in rat liver. J Pharmacokinet Biopharm. 1998;26:289-317.
PubMed 

3. Chow FS, Sharma A, Jusko WJ.  Modeling interactions between adrenal suppression and T-helper lymphocyte trafficking during multiple dosing of methylprednisolone. J Pharmacokinet Biopharm. 1999;27:559-575.
PubMed  DOI: 10.1023/A:1020974408657

4. Ferron GM, Pyszczynski NA, Jusko WJ.  Pharmacokinetic and pharmacoimmunodynamic interactions between prednisolone and sirolimus in adrenalectomized rats. J Pharmacokinet Biopharm. 1999;27:1-21.
PubMed  DOI: 10.1023/A:1020626611479

5. Ramakrishnan R, DuBois DC, Almon RR, Pyszczynski NA, Jusko WJ.  Pharmacodynamics and pharmacogenomics of methylprednisolone during 7-day infusions in rats. J Pharmacol Exp Ther. 2002;300:245-256.
PubMed  DOI: 10.1124/jpet.300.1.245

6. Ramakrishnan R, DuBois DC, Almon RR, Pyszczynski NA, Jusko WI.  Fifth-generation model for corticosteroid pharmacodynamics: application to steady-state receptor down-regulation and enzyme induction patterns during seven-day continuous infusion of methylprednisolone in rats. J Pharmacokinet Pharmacodyn. 2002;29:1-24.
PubMed  DOI: 10.1023/A:1015765201129

7. Jin JY, Almon RR, DuBois DC, Jusko WJ.  Modeling of corticosteroid pharmacogenomics in rat liver using gene microarrays. J Pharmacol Exp Ther. 2003;307:93-109.
PubMed  DOI: 10.1124/jpet.103.053256

8. Keener JP, Sneyd J. Mathematical Physiology. Berlin, Germany: Springer-Verlag; 1998. 

9. McAdams HH, Arkin A.  It’s a noisy business! Genetic regulation at the nanomolar scale. Trends Genet. 1999;15:65-69.
PubMed  DOI: 10.1016/S0168-9525(98)01659-X

10. Kepler TB, Elston TC.  Stochasticity in transcriptional regulation: origins, consequences, and mathematical representations. Biophys J. 2001;81:3116-3136.
PubMed 

11. Alm E, Arkin AP.  Biological networks. Curr Opin Struct Biol. 2003;13:193-202.
PubMed  DOI: 10.1016/S0959-440X(03)00031-9

12. Becskei A, Serrano L.  Engineering stability in gene networks by autoregulation. Nature. 2000;405:590-593.
PubMed  DOI: 10.1038/35014651

13. Ozbudak EM, Thattai M, Kurtser I, Grossman AD, van Oudenaarden A.  Regulation of noise in the expression of a single gene. Nat Genet. 2002;31:69-73.
PubMed  DOI: 10.1038/ng869

14. Cook DL, Gerber AN, Tapscott SJ.  Modeling stochastic gene expression: implications for haploinsufficiency. Proc Natl Acad Sci U S A. 1998;95:15641-15646.
PubMed  DOI: 10.1073/pnas.95.26.15641

15. Pirone JR, Elston TC.  Fluctuations in transcription factor binding can explain the graded and binary responses observed in inducible gene expression. J Theor Biol. 2004;226:111-121.
PubMed  DOI: 10.1016/j.jtbi.2003.08.008

16. Kim PM, Tidor B.  Limitations of quantitative gene regulation models: a case study. Genome Res. 2003;13:2391-2395.
PubMed  DOI: 10.1101/gr.1207003

17. Stuart RN, Branscomb EW.  Quantitative theory of in vivo lac regulation: significance of repressor packaging. J Theor Biol. 1971;31:313-329.
PubMed  DOI: 10.1016/0022-5193(71)90190-1

18. Karev GP, Wolf YI, Koonin EV.  Simple stochastic birth and death models of genome evolution: was there enough time for us to evolve? Bioinformatics. 2003;19:1889-1900.
PubMed  DOI: 10.1093/bioinformatics/btg351

19. Swain PS, Elowitz MB, Siggia ED.  Intrinsic and extrinsic contributions to stochasticity in gene expression. Proc Natl Acad Sci U S A. 2002;99:12795-12800.
PubMed  DOI: 10.1073/pnas.162041399

20. Meng TC, Somani S, Dhar P.  Modeling and simulation of biological systems with stochasticity. In Silico Biol. 2004;4:0024.
PubMed 

21. Siegel RA.  Commentary to “Linear and nonlinear system approaches in pharmacokinetics. how much do they have to offer? I. General considerations.” J Pharmacokinet Biopharm. 1988;16:667-672.
PubMed  DOI: 10.1007/BF01062017

22. Ramanathan M.  An application of Ito's lemma in population pharmacokinetics and pharmacodynamics. Pharm Res. 1999;16:584-586.
PubMed  DOI: 10.1023/A:1011910800110

23. Ramanathan M.  Pharmacokinetic variability and therapeutic drug monitoring actions at steady state. Pharm Res. 2000;17:589-592.
PubMed  DOI: 10.1023/A:1007573001055

24. Kles V, Hyrien O, Poul JM, Sanders P.  Application of pharmacokinetic, pharmacodynamic and stochastic modeling to 6-mercaptopurine micronucleus induction in mouse bone marrow erythrocytes. J Appl Toxicol. 2003;23:59-70.
PubMed  DOI: 10.1002/jat.888

25. Gillespie DT.  Exact stochastic simulation of coupled chemical reactions. J Phys Chem. 1977;81:2340-2361.
PubMed  DOI: 10.1021/j100540a008

26. Gillespie DT.  A general method for numerically stimulating the stochastic time evolution of coupled chemical reactions. J Comput Phys. 1976;22:403-434.
PubMed  DOI: 10.1016/0021-9991(76)90041-3

27. Ross IL, Browne CM, Hume DA.  Transcription of individual genes in eukaryotic cells occurs randomly and infrequently. Immunol Cell Biol. 1994;72:177-185.
PubMed 

28. Yarchuk O, Jacques N, Guillerez J, Dreyfus M.  Interdependence of translation, transcription and mRNA degradation in the lac gene. J Mol Biol. 1992;226:581-596.
PubMed  DOI: 10.1016/0022-2836(92)90617-S

29. Elowitz MB, Levine AJ, Siggia ED, Swain PS.  Stochastic gene expression in a single cell. Science. 2002;297:1183-1186.
PubMed  DOI: 10.1126/science.1070919

30. Arkin A, Ross J, McAdams HH.  Stochastic kinetic analysis of developmental pathway bifurcation in phage lambda-infected Escherichia coli cells. Genetics. 1998;149:1633-1648.
PubMed 

31. Levin MD.  Noise in gene expression as the source of non-genetic individuality in the chemotactic response of Escherichia coli.  FEBS Lett. 2003;550:135-138.
PubMed  DOI: 10.1016/S0014-5793(03)00857-3

32. Wain-Hobson S.  HIV genome variability in vivo. AIDS. 1989;3:S13-S18.
PubMed 

33. Larder BA, Darby G, Richman DD.  HIV with reduced susceptibility to zidovudine (AZT) isolated during prolonged therapy. Science. 1989;243:1731-1734.
PubMed 

34. Nijhuis M, Boucher CAB, Schipper P, Leitner T, Schuurman R, Albert J.  Stochastic processes strongly influence HIV-1evolution during suboptimal protease-inhibitor therapy. Proc Natl Acad Sci U S A. 1998;95:14441-14446.
PubMed  DOI: 10.1073/pnas.95.24.14441

35. Phillips AN, Youle M, Johnson M, Loveday C.  Use of a stochastic model to develop understanding of the impact of different patterns of antiretroviral drug use on resistance development. AIDS. 2001;15:2211-2220.
PubMed  DOI: 10.1097/00002030-200111230-00001

36. Gonze D, Halloy J, Goldbeter A.  Robustness of circadian rhythms with respect to molecular noise. Proc Natl Acad Sci U S A. 2002;99:673-678.
PubMed  DOI: 10.1073/pnas.022628299

37. Gillespie DT.  Approximate accelerated stochastic simulation of chemically reacting systems. J Phys Chem. 2001;115:1716-1733.
PubMed  DOI: 10.1063/1.1378322

38. Koyanagi S, Kuramoto Y, Nakagawa H, et al.  A molecular mechanism regulating circadian expression of vascular endothelial growth factor in tumor cells. Cancer Res. 2003;63:7277-7283.
PubMed 

39. Koyanagi S.  Optimization of the dosage schedule for sustaining intrinsic biological rhythms. Yakugaku Zasshi. 2003;123:789-797.
PubMed  DOI: 10.1248/yakushi.123.789

40. Sudo M, Sasahara K, Moriya T, Akiyama M, Hamada T, Shibata S.  Constant light housing attenuates circadian rhythms of mPer2 mRNA and mPER2 protein expression in the suprachiasmatic nucleus of mice. Neuroscience. 2003;121:493-499.
PubMed  DOI: 10.1016/S0306-4522(03)00457-3

41. Nadon R, Shoemaker J.  Statistical issues with microarray: processing and analysis. Trends Gene. 2002;18:265-271.
PubMed  DOI: 10.1016/S0168-9525(02)02665-3

42. Tian T, Burrage K. Stochastic Neural Network models for gene regulatory networks Proceedings of the 2003 IEEE Congress on Evolutionary Computation December, 2003; Canberra, Australia.