format: text

text:
# Stochastic Simulation of Gene Expression in a Single Cell

## I. Introduction

Gene expression is the central activity in a living cell. Genes are fragments of DNA molecules and provide the blueprint for the synthesis of functional molecules such as RNAs and proteins. In each cell, at any instant of time, only a subset of genes present is active in directing RNA/protein synthesis. The gene expression is “on” in such a case. There are two major steps in gene expression: transcription and translation. During transcription the sequence along one of the strands of the DNA molecule is copied or transcribed in a RNA molecule (mRNA). During translation, the sequence of the mRNA molecule is translated into the sequence of amino acids constituting a protein, i.e., a protein molecule is synthesized. Regulation of gene expression is an essential process in the living cell and determines the rates and patterns of gene expression. An in-depth understanding of gene expression and its regulation is the central focus of biology .

The processes of transcription and translation involve several biochemical reactions, the kinetics of which determine how the number of participating biomolecules changes as a function of time. In the traditional differential rate-equation approach, the time evolution of a system of chemical reactions is assumed to be continuous and deterministic. In reality, the time evolution is not a continuous process as molecular population levels in a reacting system change only by discrete integer amounts. Furthermore, the time evolution is not deterministic as the collision of molecules which brings about a chemical reaction is a probabilistic event. The deterministic rate equation approach is justified when the number of molecules of each chemical species is large compared to thermal fluctuations in the concentration. In a living cell, the number of molecules participating in different biochemical reactions is often small and there are considerable fluctuations in the reaction rates. As a result, the time evolution of the reacting system, in terms of how the number of reacting molecules changes as a function of time is stochastic rather than deterministic. There is now an increasing realization that stochasticity plays an important role in determining the outcome of biochemical processes in the cell . Stochastic effects in gene expression explain the pronounced cell-cell variation observed in isogenic populations. A cell may have the option of proceeding along one of two possible developmental pathways. The pathway selection is probabilistic and the cell fate depends on the particular choice of pathway. Thus, even a clonal population of cells can give rise to two distinct subpopulations in the course of time. The randomization of pathway choice leads to diversity and increases the likelihood of survival of organisms in widely different environments.

The issue of stochasticity (randomness or noise) and its effect on cellular processes as well as on the operation of synthetic devices like genetic switches, has been addressed in several theoretical studies . A complete understanding of cellular processes requires an appreciation of events at the level of an individual cell and subsequent extrapolation to an ensemble of cells. Recent experimental advances have made it possible to study processes within a single cell unmasked by ensemble averaging . The simplest event one can study at the single cell level is that of the expression of a reporter gene such as *lacZ* and *GFP*. In the former case, the end product is an enzyme $`\beta `$-galactosidase which is capable of hydrolyzing a noncolored substrate to a colored product. In the latter case, the protein itself is fluorescent. Hence, the gene expression can be directly studied either colorimetrically or fluorometrically at the level of an individual cell. Recent experiments using such techniques, provide evidence that gene expression occurs in abrupt stochastic bursts at the level of an individual cell . Two very recent experiments  provide direct evidence of stochasticity in gene expression. In both the experiments, a quantitative measure of the noise associated with gene expression has been obtained. The noise has both intrinsic and extrinsic components. Intrinsic noise is the difference in protein synthesis which arises when two identical copies of a gene are expressed under the same conditions. Extrinsic noise occurs due to fluctuations in the cellular components required for gene expression. The experiments provide a quantitative framework for the characterization of noise in gene regulatory networks.

Some earlier experiments on both prokaryotic and eukaryotic cells have provided evidence of the so-called “all or none” phenomenon in gene expression . This implies that in an individual cell, gene expression is either low/off or has a high value. In an ensemble of cells the protein levels are distributed in a bimodal manner, a large fraction of cells synthesize proteins at a low (may be zero) level or produce them at a high level. Most of the experiments require the presence of inducers/enhancers to observe bimodality. There is strong experimental evidence that inducer/enhancers increase the number of expressing cells but not the level of expression per cell . The process of gene expression is analogous to a binary switch which can be in “on” and “off” positions. Inducer/enhancer molecules make it favorable for the switch to be in the “on” position. Some theories have been proposed so far to explain the so-called “all or none” phenomenon in prokaryotic gene expression. The theories are mostly based on an autocatalytic feedback mechanism , synthesis of the gene product gives rise to the transport or production of inducer molecules which in turn promote further gene expression. In section II of this paper, we propose a model of gene expression and show that in the presence of a sufficient number of inducer molecules in a cell, gene expression in that cell occurs at a high level. In cells where inducer molecules are absent or are few in number, gene expression occurs at practically zero level. The role of autocatalysis in the “all or none” phenomenon is also commented upon. The method employed for the study is that of stochastic simulation based on the Gillespie Algorithm (GA) . The GA provides a stochastic realization of the temporal pattern of gene expression and is more realistic and accurate than the deterministic differential rate equation approach. Our model of gene expression includes the major biochemical reactions involved in transcription and translation. In section III, we explore the parameter space of the model and obtain different temporal patterns of gene expression. One parameter region of particular interest corresponds to stochastic flips between the states 1 and 2 at random time intervals. In state 1, the protein level is zero. This is an example of a binary switch which makes stochastic transitions between the states 1 and 2 and the temporal process is analogous to a two-state jump phenomenon. The effect of changing the reaction parameters on the temporal patterns of gene expression is further studied. Some of the results can be understood in the framework of a simple mathematical model.

## II. Stochastic model of gene expression

We consider a single gene. The gene is transcribed into mRNA by an enzyme called RNA polymerase (RNAP). The process is initiated with the binding of RNAP to a site called promoter, usually near the beginning of the transcribed sequence. Expression of most genes are regulated at the level of transcription and more specifically during the initiation of transcription, that is, before the first phosphodiester bond is formed. Regulation of transcription initiation is achieved by the binding of a regulatory protein (R) to an overlapping segment of DNA (called operator O) resulting in a turning off of mRNA production. RNAP and the regulatory R molecules are mutually exclusive. If RNAP binds to the promoter region first, it prevents the binding of R to the operator region and vice versa. An inducer molecule (I) may bind to R both when R is free and when R is bound to the operator O. In the later case, the complex of I and R detaches from the operator. As long as R is forming a bound complex with I, it is unable to bind at O and so cannot function as a regulatory protein. R regains its activity when the inducer dissociates from the I\_R complex and R is able to occupy the operator region once more. The biochemical reactions considered in the model of gene expression are:
Reaction 1:

$$\text{O}+\text{R}=\text{O}_{}\text{R}$$
(1)
Regulatory molecule R binds to the operator region O to form the bound complex $`\text{O}_{}\text{R}`$.
Reaction 2:

$$\text{O}_{}\text{R}\text{O}+\text{R}$$
(2)
Bound complex O\_R dissociates into free R and O.
Reaction 3:

$$\text{P}+\text{RNAP}=\text{P}_{}\text{RNAP}_{cc}$$
(3)
RNA polymerase (RNAP) binds to the promoter region P forming the closed complex P\_RNAPcc.
Reaction 4:

$$\text{P}_{}\text{RNAP}_{cc}\text{P}+\text{RNAP}$$
(4)
The closed complex dissociates into free RNAP and P.
Reaction 5:

$$\text{P}_{}\text{RNAP}_{cc}\text{P}_{}\text{RNAP}_{oc}$$
(5)
Isomerization of closed to open complex P\_RNAPoc occurs. The open complex is the activated form of the RNAP-promoter complex.
Reaction 6:

$$\text{P}_{}\text{RNAP}_{oc}\text{TrRNAP}+\text{RBS}+\text{P}$$
(6)
RNAP clears the promoter region and synthesis of the mRNA chain starts. The appearance of the ribosome binding site RBS occurs at the beginning of the mRNA chain. TrRNAP denotes transcribing RNA polymerase.
Reaction 7:

$$\text{TrRNAP}\text{RNAP}$$
(7)
RNAP completes transcription and is released from DNA.
Reaction 8:

$$\text{RBS}+\text{Ribosome}\text{RibRBS}$$
(8)
Ribosome binds to RBS and RibRBS denotes the bound complex.
Reaction 9:

$$\text{RibRBS}\text{RBS}+\text{Ribosome}$$
(9)
Ribosome dissociates from the bound complex RibRBS.
Reaction 10:

$$\text{RBS}\text{degradation}$$
(10)
RBS degrades due to the binding of RNAseE at RBS. This binding event is not considered separately.
Reaction 11:

$$\text{RibRBS}\text{EIRib}+\text{RBS}$$
(11)
RBS is cleared and the ribosome EIRib initiates translation of mRNA chain.
Reaction 12:

$$\text{EIRib}\text{protein}$$
(12)
Protein synthesis by transcribing ribosome is completed.
Reaction 13:

$$\text{protein}\text{degradation}$$
(13)
Degradation of protein product occurs.
Reaction 14:

$$\text{I}+\text{R}\text{I}_{}\text{R}$$
(14)
Inducer molecule binds to free regulatory molecule R. I\_R is the bound complex of I and R.
Reaction 15:

$$\text{I}_{}\text{R}\text{I}+\text{R}$$
(15)
Bound complex dissociates into free I and R.
Reaction 16:

$$\text{O}_{}\text{R}+\text{I}\text{I}_{}\text{R}+\text{O}$$
(16)
Inducer binds to bound complex $`\text{O}_{}\text{R}`$, the complex $`\text{I}_{}\text{R}`$ detaches and the operator region O is freed. We emphasize that for many of the steps described above, alternate mechanisms exist. However, the mechanism described here is consistent with many gene regulatory systems.

Reaction schemes (1) - (13) are based on those considered in Refs.. Transcription and translation are tightly coupled in prokaryotes. As soon as RNAP leaves the promoter region, the 5 end of the mRNA chain, containing the RBS, is available for ribosome binding (Reaction 6). This implies that the mRNA chain need not be completely synthesized to allow for ribosome binding (Reaction 8), protein synthesis by translating ribosome (Reaction 12) and mRNA degradation at RBS (Reaction 10). Following Ref., the number of mRNA molecules at any instant of time is given by the sum of the numbers of RBS and RibRBS, the bound complex of ribosome and RBS. The functional degradation of mRNA starts at the moment RNAseE binds to the RBS.

We now give a brief description of the Gillespie Algorithm (GA) . Suppose there are N chemical species participating in M chemical reactions. Let X(i), i = 1, 2, 3 . . . . , N denote the number of molecules of the i th chemical species. Given the values of X(i), i = 1, 2, 3, . . . , N at a time $`t`$, the GA is designed to answer two questions: (1) when will the next reaction occur? and (2) what type of reaction will it be? Let the next reaction occur at time $`t`$$`+`$$`\tau `$. Knowing the type of reaction, one can adjust the numbers of participating molecules in accordance with the reaction schemes. Thus, with repeated applications of the GA, one can keep track of how the numbers X(i)’s change as a function of time due to the occurrence of M different types of chemical reactions. Each reaction $`\mu `$ ($`\mu `$ = 1, 2, . . . , M) has a stochastic rate constant $`c`$μ associated with it. This rate constant has the following interpretation:

$`c`$μ$`dt`$ = probability that a particular combination of reactant molecules participates in the $`\mu `$ th reaction in the infinitesimal time interval ($`t`$, $`t+dt`$). Let $`h`$μ be the number of distinct molecular combinations for the $`\mu `$ th reaction. Then

$`a`$μ$`dt`$$`=`$$`h`$μ$`c`$μ$`dt`$ = probability that the $`\mu `$ th reaction occurs in the infinitesimal time interval ($`t`$, $`t+dt`$).
Let P($`\tau `$,$`\mu `$) $`d`$$`\tau `$ be the probability that the next reaction is of type $`\mu `$ and occurs in the time interval ($`t`$\+ $`\tau `$, $`t`$+$`\tau `$ \+ d$`\tau `$). It is straightforward to show that

$$\text{P}(\tau ,\mu )=a_\mu \mathrm{exp}(a_0\tau )$$
(17)
where

$$a_0=\underset{\nu =1}{\overset{\text{M}}{}}a_\nu \underset{\nu =1}{\overset{\text{M}}{}}h_\nu c_\nu $$
(18)
What is needed now is to generate a random pair ($`\tau `$,$`\mu `$) according to the probability distribution (17). Let r1 and r2 be two random numbers obtained by invoking the standard unit interval uniform random number generator. One can then show that $`\tau `$ and $`\mu `$ are obtained as

$$\tau =(\frac{1}{a_0})\mathrm{ln}(\frac{1}{r_1})$$
(19)
and $`\mu `$ is taken to be the integer for which

$$\underset{\nu =1}{\overset{\mu 1}{}}a_\nu <r_2a_0\underset{\nu =1}{\overset{\mu }{}}a_\nu $$
(20)
A rigorous proof of the formulae (19) - (20) is given in Ref.. Once $`\tau `$ and $`\mu `$ are known, the time evolution of the reacting system is specified. The stochastic rate constant $`c`$μ is related to the more familiar deterministic reaction rate constant $`k`$μ through simple relations. In the case of first order reactions, both constants have the same value. In the case of second order reactions, the rate constant is divided by the volume of the system. We have applied the GA to our model of gene expression involving N = 15 types of biomolecules participating in M = 16 biochemical reactions. In the initial state (time $`t`$=0), the number of free operator and promoter sites is 1. The number of R, I, RNAP, and ribosome molecules is $`N_R=`$ 20, $`N_I=`$ 20, $`N_{RNAP}=`$ 400 and $`N_{Rib}=`$ 350 respectively. The number of all the other biomolecules is set to zero at $`t`$ = 0. The simulation time is up to 2000s, i.e., less than the cell generation time typically in the range 2000-3000s. As already mentioned, knowledge of $`\tau `$ and $`\mu `$ enables one to calculate the numbers of biomolecules at time $`t`$+$`\tau `$ and in this way, through repeated applications of the GA, one can keep track of how the different numbers change as a function of time. Fig.1 shows the number of protein molecules present in the system as a function of time. The stochastic rate constants of the sixteen reactions are c1 = 0.5, c2 = 0.004, c3 = 0.02, c4 = 0.001, c5 = 0.8, c6 = 0.9, c7 = 0.08, c8 = 0.01, c9 = 0.001, c10 = 0.3, c11 = 0.8, c12 = 0.7, c13 = 0.003, c14 = 0.7, c15 = 0.001 and c16 = 0.3 respectively. The topmost curve corresponding to $`N_I`$ shows that the protein number reaches a steady level with fluctuations around the mean. The dashed curve at the bottom of the figure shows the number of proteins as a function of time when the number of inducer molecules $`N_I`$ (=3) is much less than $`N_R`$. The protein number in this case is very small, practically zero. The other curves corresponds to values of $`N_I`$ closer to $`N_R.`$ As $`N_I`$ approaches $`N_R`$, protein levels are intermediate between high and low. Figure 2 shows the number of protein molecules as a function of time in a different parameter region.

The “all or none” phenomenon has been observed in both prokaryotic and eukaryotic cells. In prokaryotes, genes are often arranged in operons,i.e, sets of contiguous genes which include structural and regulatory sequences. A well-known example is that of the *E. coli* lactose (*lac*) operon . Lac operon consists of three structural genes (z, y and $`a`$) which code for the three proteins: $`\beta `$-galactosidase, the enzyme that catalyzes the hydrolysis of lactose to glucose and galactose; permease, a carrier protein responsible for membrane transport of lactose into the cell and a third protein transacetylase. The lac operon contains three regulatory sequences, i, P and O, which control the transcription of mRNA leading to the synthesis of the three proteins. The sequence i corresponds to a gene lacI which is transcribed continuously to synthesize a repressor protein, lac repressor, at a low level. Lac repressor binds to the operator sequence O and prevents the transcription of the genes z, y and $`a`$ so that the $`\beta `$-galactosidase enzyme and the permease molecules are not produced. If the bacterium is to grow on lactose (milk sugar) which acts as its carbon source, $`\beta `$-galactosidase must be made available to split the sugar into glucose and galactose. The breakdown product of lactose act as an inducer molecule. The inducer attaches to the repressor molecule, causing it to release the DNA so that transcription of the structural gene is possible. The repressor is freed of the inducer when the lactose supply is exhausted and switches off the expression of the structural genes once more. Our simple model of gene expression incorporates some of the key features of the lac operon. The role of the lac repressor is played by the regulatory molecule R though its synthesis is not explicitly considered. There is a single gene in our model analogous to the structural gene z expressing the enzyme $`\beta `$-galactosidase. The inducer molecule I acts in the same manner as in the case of the lac operon.

The experimental observation of the “all-or-none” phenomenon in the lac operon has been attributed to autocatalytic feedback mechanisms . At low inducer concentrations, some of the bacterial cells synthesize protein at the full rate whereas the other cells are in the “off” state. When inducer is added to the colony of bacterial cells, simultaneous production of the $`\beta `$-galactosidase enzyme and permease molecules occurs. The permease molecules transport lactose into the cell raising the internal inducer concentration which in turn promotes the production of more $`\beta `$-galactosidase and permease. Thus an autocatalytic feedback process is at work and within a short time after the appearance of the first permease molecules, the bacterial cell becomes fully induced synthesizing the structural proteins at maximum rate. Siegele and Hu  carried out experiments on gene expression from plasmids containing the araBAD promoter in the presence of subsaturating concentrations of the inducer arabinose. Again, as in the case of the lac operon, it has been suggested that an autocatalytic induction mechanism, due to the accumulation of inducer molecules by transport, is at work. However in all the experiments, full induction of cells has been observed even in the absence of autocatalysis, i.e., when the inducer availability is not linked to that of synthesized protein molecules like permease. Our model of gene expression does not include an autocatalytic feedback process and a detailed analysis of the simulation results shows that if the number of inducer $`(`$$`\text{I})`$ molecules, $`N_I`$, is greater than or equal to the number of regulatory $`(\text{R})`$ molecules, $`N_R`$, in a cell, the cell reaches a steady state which is state 2 (high protein level). If $`N_I`$= $`0`$ or $``$ $`N_R`$, the cell is in state 1 (low/zero level) but as $`N_I`$ approaches $`N_R`$ , protein levels intermediate between high and low are obtained. The “all or none” phenomenon becomes more pronounced in the presence of the autocatalytic induction (AI) mechanism. If $`N_I`$ is originally small in a cell, the AI mechanism leads to a rapid increase in $`N_I`$ . When $`N_I`$ is $``$ $`N_R`$, the cell is in state 2 in the steady state. The magnitude of the protein level in state 2 is independent of the value of $`N_I`$ in agreement with experimental results. Autocatalysis is responsible for a sharp bimodal distribution in protein levels. In the absence of autocatalysis, the magnitude of the protein level in state 2 is the same as in the case of autocatalysis but the bimodal distribution become less sharp due to intermediate protein levels in a fraction of cells. Our simulation results are in agreement with experimental observations . Experiments  on single mammalian cells (eukaryotic cells) have provided further evidence of bimodality in the distribution of protein levels in an ensemble of cells. Again, one finds that the amount of enhancer affects the number of expressing cells but not the level of expression. In other words, the enhancer increases the probability rather than the rate of transcription. In the experiments carried out by Zlokarnik et al. , the reporter gene synthesizes the protein $`\beta `$-galactosidase. In unstimulared cells, the number of these proteins is low, in the range 150 - 300. Under the action of the stimulating agent carbachol, rapid conversion to a state, in which 15000 - 20000 $`\beta `$-lactonase molecules are present in a cell, is obtained. The major conclusion of the single cell experiments mentioned above is that in the systems considered, the cellular state is bistable. A cell can exist in two stable steady states: gene expression “off/low” and gene expression “on” with a high level of protein production. Addition of inducer/enhancer to the system increases the fraction of cells in the “high” state. In the case of eukaryotic systems, however, the mechanism of enhancer action is not well understood. Enhancers have been suggested to give rise to two major types of response. Enhancers increase the rate of transcription through mainly enhancing the rate of close to open complex formation of RNAP bound to the promoter region . The second type of response is of the “all or none” type . Enhancers in this case do not increase the rate of transcription but increase the fraction of cells in the high state (state 2). The “all or none” phenomenon observed in some eukaryotic systems  dose not involve autocatalysis explicitly. Thus, a more general mechanism than in the case of prokaryotic systems is required to explain the bimodal distributions in protein levels in an ensemble of cells.

## III. Patterns of gene expression

We now consider the model of gene expression in the absence of inducer molecules. The values of M, the total number of reactions and N, the number of different types of biomolecules, are both thirteen. To explore the full parameter space, one has to treat the thirteen stochastic rate constants corresponding to the same number of reactions as variables. Experimental results wherever they are available, show that the usual rate constants, $`k`$μ’s, to which the stochastic rate constants $`c_\mu `$’s are related, can vary over a wide range depending on the type of gene and the nature of the cellular environment . Since the exploration of the full parameter space is a daunting task, we report on some of the more general patterns of gene expression in a restricted subspace. The effect of changing the stochastic rate constants on specific patterns is also studied. We have included an optional feature in our model of gene expression, namely, cooperative binding of RNAP to the promoter region P. The possibility of such a binding has been suggested earlier in a simplified probabilistic model of gene expression . In the present model, cooperative binding implies that the rate constant for the binding of RNAP at P is enhanced by a factor q if the binding event is immediately preceded by the Reaction 6 in which a RNAP clears the promoter region. Cooperative binding of proteins to DNA is now well established. In most cases of regulatory proteins, the binding cooperativity is mediated through protein-protein interaction although increasing evidence of DNA mediated effects have been reported . Some recent experiments on the prokaryotic system E. coli have shown that the transcriptional activity of the promoter is intrinsically sensitive to the superhelical density of the DNA template . In fact as pointed out by McClure, supercoiling gives rise to considerably more diversity in the patterns of promoter strength (the ability to bind weakly or strongly) than do mutations of auxiliary proteins. There is now experimental evidence that transcription generates increased negative supercoiling through several hundred base pairs. In principle, this can facilitate the binding of RNAP to the DNA or decrease the energy of activation required for the isomerization of RNAP-promoter complex from closed to open form. Thus, it is entirely plausible and likely that active transcription downstream of the promoter site may lead to increased binding of RNAP (cooperative binding) and enhanced rate of open complex formation (stochastic rate constants $`c_5`$).

We now describe three major types of gene expression pattern as a function of time. As before, the GA is applied to the system of thirteen biochemical reactions constituting the processes of transcription, translation and the regulation thereof. In the gene expression pattern designated as Type A, the protein production occurs in abrupt stochastic bursts. A variable number of proteins is produced in each burst. The type A pattern of gene expression has been observed experimentally and has been attributed to stochastic effects. Figures 3 and 4 show two such patterns of protein production for $`q=1`$ (no cooperative binding).

In the Type B pattern, the protein level reaches a steady state with fluctuations around the mean (Figures 5 and 6). This type pattern is quite common and is routinely observed in experiments.

The Type C pattern of gene expression has an interesting structure. As in the case of the Type A pattern, protein production occurs in stochastic bursts, i.e., at random time intervals. The bursts may be of various durations but in each burst, the protein number attains the same level (with attendant fluctuations) in a very short time. Similarly, the decay of the protein level from high to zero occurs in a small time interval. Figures 7(a) and 7(b) show the patterns of mRNA and protein production as a function of time. The number of regulatory (R) molecules ($`N_R`$), RNAP ($`N_{RNAP}`$) and ribosome ($`N_{Rib}`$) is 10, 400 and 200 respectively.

We now discuss the physical origin of the Type C pattern. The stochastic rate constant c3 for RNAP binding (Reaction 3) is considerably higher than that of the binding of R at the operator site O (Reaction 1). The initial number of RNAP molecules is also larger than that of R molecules. As explained in section II, $`a_\mu `$$`dt`$ $`=`$ $`h_\mu c_\mu dt`$ is the probability that the $`\mu `$ th reaction occurs in the infinitesimal time interval ($`t`$ , $`t+dt`$ ). The number of distinct molecular combinations $`h_\mu `$ for the $`\mu `$ th reaction is 10 and 400 for the Reactions 1 and 3 respectively. The corresponding stochastic rate constants have the values c1= 0.01 and c3 = 0.7. Thus, Reaction 3 is more probable than Reaction 1. Note that the stochastic rate constants c5 , c6 are considerably high. Reactions 5 and 6 are associated with the transcription process, namely, isomerization of the closed complex of RNAP bound to the promoter region P to the open complex and subsequent clearance of the promoter region by RNAP. After the binding of a RNAP to P, a host of factors including high values of some of the relevant rate constants, leads to a sharp rise in the number of proteins to a level determined by the transcription, translation and protein degradation rates. The protein level is maintained over a time interval due to the balancing of the rates of synthesis and degradation. Binding of the R molecule to O, though less probable than that of RNAP at P, can occur with a finite probability. Once the R molecule is bound to O, it continues to remain bound for some time as the dissociation rate (c2 = 0.004) is low. This prevents the binding of a RNAP to P during the time interval in which R stays bound to O, leading to a sharp fall in the number of proteins to zero. The subsequent dissociation of the R molecule from the operator O, followed by the binding of a RNAP to the promoter region P, tilts the balance in favour of state 2. The cellular state thus flips from state 1 (no. of proteins zero) to state 2 (no. of proteins high) and vice versa at random time intervals dictated by stochastic binding and dissociation events at O. In section II, we have discussed a bimodal distribution in protein levels due to the “all or none” phenomenon in an ensemble of cells. In the present case, a bimodal distribution in protein levels is obtained in an ensemble of cells if protein levels are measured at particular instant of time. The heights and widths of the two peaks may change as a function of time but the distribution remains bimodal.

Figure 8 shows the distribution of the number N(m) of cells expressing a fraction m of the average number of proteins at a particular instant of time. The total number of cells is 1000. The values of the stochastic rate constants are the same as in Figs. 7(a) and 7(b). Type C pattern of gene expression is also obtained for lower values of $`c_5`$ and $`c_6`$ (Fig. 9) though better quality patterns are obtained for high values of the rate constants. Figures 10(a) and 10(b) show the temporal variations of the mRNA and protein numbers with stochastic rate constants the same as in Figs. 7(a) and 7(b) but the enhancement factor q, associated with cooperative RNAP binding, has been raised from 1 to 10. Comparing the two sets of Figures, one concludes that cooperativity increases the duration of state 2 (“high” level). The magnitude of the mean level remains unchanged. In an ensemble of cells, a greater fraction of cells is in state 2 than in the earlier case. The origin and nature of bimodal distribution in protein levels are different for the model system considered in section II and the type C pattern of gene expression. In the latter case, no autocatalytic feedback is necessary to obtain a bimodal distribution. Positive (autocatalytic) feedback mechanism has been invoked to explain the “all-or-none” phenomenon in prokaryotic and eukaryoticsystems. Experimental reports of the phenomenon in some eukaryotic systems suggest that autocatalytic (positive) feedback is not essential for the obsevance of the phenomenon. In these systems, activation to the high protein level is enhancer mediated. The origin of bimodal distribution in protein levels is, however, yet to be elucidated.

Figure 11 shows the number of proteins (solid line) and mRNA molecules (dotted line) as a function of time to make a simultaneous comparison of the rise and decay of protein and mRNA levels.Note that after the number of mRNA molecules become zero, there is a time delay before the number of proteins falls to zero. Even when there are no mRNA molecules in the system, some proteins remain which degrade to zero level at a decay rate lower than that of the mRNA molecules. In Fig. 12, a pattern of gene expression is shown in which the protein number never falls to zero and a variable number of proteins is synthesized as a function of time.

We have further checked how robust pattern C ( the gene expression pattern shown in Figure 7) is when the different stochastic parameters are changed.

We change one parameter at a time keeping all the other parameter values the same as in Fig. 7. As c1 decreases from 0.01 (c1 = 0.01 in Fig. 7), the duration of state 2 increases and ultimately state 2 becomes the steady state (Figs. 13(a) and 13(b)), i.e., a Type B pattern is obtained. As c1 increases from 0.01, the Type C pattern is still obtained (Fig. 13(c), c1 = 0.02) but for higher c1 values, say, c1 = 0.03, the pattern of gene expression becomes of Type A (Fig. 13(d)). If c2 is varied, then as c2 decreases from 0.004 (Fig. 7), the total duration of state 2, for $`t`$ in the range 0 - 2000s, decreases (Fig. 13(e) with c2 = 0.008). If c4 is varied, one finds that the Type C pattern of gene expression is obtained over a wide range of values. The greater the values of c5 and c6 , the higher is the level of proteins attained in state 2. The best Type C patterns are obtained if the values of c5 and c6 are high. Figs. 13(f) and 13(g) correspond to c5= 2, c6 = 20 and c5= 20, c6 = 2 respectively. If c7 is varied, the Type C pattern is obtained over a wide range of values. If c8 increases from the value 0.01 (Fig. 7), Type C pattern becomes Type A (Fig. 13(h), c8 = 0.015). If c8 decreases from the value 0.01, the total duration of state 2 increases(Fig. 13(i)). If c10 decreases from the value 0.2 (Fig. 7), the Type C pattern is lost (Fig. 13(j), c10 = 0.1). If c10 increases from 0.2 the Type C pattern is obtained over a wide range of values of c10. The best type of Type C pattern is obtained for similar values of c11 and c12.

We have also studied the effect of changing the number $`N_R`$, $`N_{RNAP}`$ and $`N_{Rib}`$ of R, RNAP and ribosome molecules respectively. As $`N_R`$ increases from 10 (Fig. 7), the Type C pattern is gradually lost. Reducing $`N_R`$ increases the total duration of state 2. If$`N_R`$= 0, state 2 becomes the steady state. If $`N_{RNAP}`$ is reduced from 400 (Fig. 7) the Type C pattern worsens gradually. If $`N_{RNAP}`$ increases beyond 400, the Type C pattern is still obtained. The magnitude of the protein level in state 2 remains more or less the same. If $`N_{Rib}`$ decreases from 200, the total duration of state 2 increases.

Some general conclusions that can be made, on the basis of the stochastic simulation of gene expression patterns, are as follows. The protein level attained in state 2 is not affected by changes in the rate constants c1, c2, c4, c7, c11, c12 and the enhancement factor q. The rate constants c5, c6, c8, c10 and c13 determine the magnitude of the protein level. If c1 decreases and c2 increases, the total duration T of state 2 in the time interval 0-2000s increases. Transition from Type C $``$Type B can occur by decreasing c1(Fig. 13(b)), increasing q and by making the values of c5 and c6 unequal (Fig. 13(g)). Transition from Type C $``$ Type A pattern can occur by increasing c1 (Fig. 13(d)), c8 (Fig. 13(h)) or by decreasing c10 (Fig.13(j)). Type C patterns are favourable for high values of c5 and c6 with c5$``$ c6 and also c11 $``$ c12.

Kepler and Elston have proposed a simple mathematical model of gene expression with no feedback. In their model, if an activator molecule occupies the operator region, protein production occurs at a rate $`\alpha _1.`$ If the operator is unoccupied by the activator molecule, protein production occurs at a lower rate $`\alpha _0.`$ The model does not include the intermediate steps of protein synthesis like transcription, ribosome binding for the initiation of translation etc. The model, however, provides lots of physical insight on stochastic effects in the form of fluctuations in the discrete states of the operator (unoccupied/occupied) on gene expression. Chemical reactions that change the state of the operator are termed operator fluctuations.Approximations to the dynamics were made for the cases in which the protein number is large or the operator fluctuations are fast. In the first case, the effective rate of protein synthesis fluctuates randomly in time between high (synthesis rate $`\alpha _1`$) and low (synthesis rate $`\alpha _0`$) levels. In the latter case, the fluctuations are effectively averaged out over larger time scales. In our model, the major biochemical reactions/events in protein synthesis have been included and instead of activators we have regulatory molecules which act as repressors. The simulation based on GA provides an accurate knowledge of the microscopic origins of the different types of temporal gene expression patterns. A simple mathematical model provides insight on the origin of different patterns and transition from one type pattern to another. Let $`m`$ denote the number of proteins at time t devided by the maximum number of proteins. The equation

$$\frac{dm}{dt}=xm$$
(21)
describes the rate of change in the number of proteins. The possible values of $`x`$ are $`0`$ and $`1`$ so that in the steady state $`m`$ can be either $`1`$ (high level, state 2) or $`0`$ (low level, state 1). The variable $`x`$ randomly switches between the two states and the transition rate from state 1 $``$ state 2 is $`r_1`$ and that from state 2 $``$ state 1 is $`r_2`$. Let $`P_j(m=M,t)`$, $`j=1,2`$, be the probabilities of being in the state $`j`$. One can then write down the master equations

$$\frac{P_0}{t}=\frac{(MP_0)}{t}r_1P_0+r_2P_1$$
(22)
$$\frac{P_1}{t}=\frac{((1M)P_1)}{t}r_2P_1+r_1P_0$$
(23)
The steady state distribution of $`P(m,t)=P_0(m,t)+P_1(m,t)`$ is given by

$$P(m=M)=AM^{r_11}(1M)^{r_21}$$
(24)
where $`A`$ is the normalization constant. Figs. 14(a), (b) and (c) show the distribution $`P(m)`$ versus $`m`$ for different values of $`r_1`$ and $`r_2`$. The distribution in Fig. 14(a) corresponds to the Type C pattern of gene expression (the magnitude of the high protein level in state 2 is normalized to the value $`1`$) in the steady state and is bimodal in nature. Fig. 14(b) describes the Type B pattern of gene expression and there is a single peak corresponding to the steady state level, the magnitude of which is normalized to $`1`$. Fig. 14(c) shows a broad distribution in protein levels. Fig. 12 shows a pattern of gene expression based on simulation of the detailed model which gives rise to a probability distribution of protein levels similar to that shown in Fig. 14(c). The transition rates $`r_1`$ and $`r_2`$ are functions of the different stochastic rate constants though the actual functional relationship is yet to be worked out. The rate $`r_1`$ depends dominantly on $`c_2,`$ $`c_5,`$ $`c_6,`$ $`c_8`$ and $`c_{10}`$ whereas the rate $`r_2`$ is mainly determined by $`c_1`$ and $`c_{13}.`$ The simple model illustrates how the different probability distributions arise and the transition from one type to another occurs as the transition rates $`r_1`$ and $`r_2`$ are changed. For slow transition rates, the states 1 and 2 can be distinguished and the probability distribution $`P(m)`$ is bimodal. For fast transition rates, $`m`$ has values intermediate between $`0`$ and $`1`$ and the bimodality is smeared into a single broad distribution. For large $`r_1`$ and small $`r_2`$, the peak is around the high value $`1`$. The transitions from one type of probability distribution to another as $`r_1`$ and $`r_2`$ are varied are consistent with the changes in the nature of the temporal patterns of gene expression (Figs.12 and 13), brought about by changes in the various stochastic rate constants on which the transition rates $`r_1`$ and $`r_2`$ of the simpler model depend.

## IV. Concluding remarks

In this paper, we have considered two models of gene expression in a single cell. In the first model, the effect of inducer molecules is explicitly considered and the cellular kinetics are described by a set of sixteen biochemical reactions. In the second model, the cellular kinetics are described by a set of thirteen biochemical reactions. Both the models include the basic steps of transcription and translation as well as transcriptional regulation through the binding of a regulatory molecule R to the operator region, with R acting as a repressor of transcription. Additionally, in the second model a transcriptionally generated cooperativity factor q in the binding of RNAP has been introduced as an optional feature. The patterns of gene expression as a function of time have been determined in both the models with the highly accurate stochastic simulation method based on the Gillespie Algorithm. Analytic approaches  are possible in studying the temporal evolution but in these cases the detailed biochemical reactions are lumped together into a few effective processes. This makes it possible to determine the temporal evolution of the system in a chemical Master Equation approach which is stochastic in nature. Alternatively, the effect of stochasticity can be taken into account by the inclusion of “noise” terms in the differential rate equations. The GA provides a more detailed picture of the kinetics though the computational efforts required in its implementation are considerable. In section II, we have studied gene expression in the presence of inducer molecules. The major result of this study is that when the number of inducer molecules, $`N_I`$, is greater than or equal to the number of regulatory molecules, $`N_R`$, the cellular steady state is state 2 in which the proteins are synthesized at a high level. When $`N_IN_R`$ , the cellular state is state 1 in which the protein level is low/zero. As $`N_I`$ approaches $`N_R`$, protein production occurs at intermediate levels. Fig.1 provides an example of how drastic the effect of even a single molecule can be on gene expression. There is a considerable difference of protein levels attained in the steady state when $`N_I`$ is changed from $`19`$ to $`20`$. This single molecule effect can be observed over a wide region of parameter space and may give rise to threshold phenomena. These results are new and have not been reported earlier. McAdams and Arkin have pointed out that even a single molecule can switch the biochemical state of a cell. Togashi and Kaneko  have studied an autocatalytic reaction system with a small number of molecules and shown that due to the nonlinear dynamics, amplification of small changes can give rise to single molecule switches.

The results of our simulation explain why an autocatalytic induction mechanism gives rise to the “all-or-none” phenomenon observed in cells in the presence of inducer molecules. Due to autocatalytic induction, the small number of inducer molecules initially present in a cell is quickly amplified so that the condition $`N_IN_R`$ is satisfied in a short interval of time. The cell then exists in state 2 with a high protein level. Other cells in which $`N_I`$ is zero exist in state 1 with zero protein level. In the presence of a subsaturating concentration of inducer molecules, the distribution of protein levels in an ensemble of cells is bimodal. In the absence of autocatalytic induction, protein production occurs at low, high as well as at intermediate levels so that the bimodal distribution gets smeared. This is in keeping with experimental observations.

The model studied in section III does not include inducer molecules. Simulation based on the GA shows the existence of three types of pattern of gene expression as a function of time. In the Type A pattern, protein synthesis occurs in abrupt stochastic bursts and a variable number of protein is produced in each burst. There is considerable experimental evidence  for this type of gene expression and in section III, some examples of this type of pattern have been given (Figs.3 and 4). In Type B pattern, protein levels reach a steady state (Figs. 5 and 6). This is very common type of pattern observed in gene expression experiments.

In the Type C pattern (Fig. 7), the cellular state makes random transitions between states 1 and 2 . This is a different manifestation of the “all or none” phenomenon and measurement levels in an ensemble of cells is bimodal (Fig. 8). Type C patterns of gene expression have been obtained by Kepler and Elston  in their study of model gene expression systems (see Fig. 4 of Ref. ) using the Master Equation Approach. In these models detailed biochemical reactions have been replaced by a few effective processes so that the models are mathematically tractable. Kepler et al. have specifically considered the effect of fluctuations in the state of the operator on gene expression. Their conclusion is that the operator fluctuations can induce bistability in parameter regions which gives rise to monostability in the deterministic, i.e., the zero noise limit or destroy bistability if it exists in the noise-free case. In the deterministically bistable region, the gene acts like a genetic switch and external noise/perturbation is needed to flip the switch from one state to the other. This is the principle of operation behind the noise-based switches and amplifiers for gene expression proposed by Hasty et al . If a system is stochstically bistable, the fluctuations in the system flip the switch between the two states (say, state 1 and state 2) at random time intervals. As in Ref. , the operator fluctuations have been explicitly considered in the models studied in Sections II and III and it has been shown that in certain parameter regimes, stochastically bistable behaviour corresponding to Type C pattern of gene expression is obtained. Our stochastic simulation results, obtained by using the highly accurate Gillespie Algorithm, provide a microscopic basis as well as quantitative estimates of the different rate constants for obtaining Type C pattern of gene expression. The simulation method can be used to study a large number of reactions which is not feasible in the formalism of a mathematical model. The significant omission in Ref.  is that no distinction has been made between transcription and translation. The models describe direct translation from gene into protein. As pointed out in the paper, the simplification may have considerable impact on cellular phenomena. For example, the synthetic repressilator network  would not oscillate if a time delay between transcription and translation were absent. In our model of gene expression, all the major biochemical reactions involved in transcription and translation have been explicitly taken into account and the time delay is clearly seen in Fig. 11. There is some experimental evidence  of transcriptionally generated cooperative binding of RNAP to the promoter region. For non-zero valus of the cooperativity factor $`q`$ in our model the duration of state two (high protein level) is found to increase. The result is new and the effect of cooperative RNAP binding on gene expression needs to be investigated in greater detail.

We have further studied the effect of changing the various stochastic rate constants on the temporal pattern of gene expression. Transitions involving Type C $``$ Type B and Type C $``$ Type A patterns of gene expression have been obtained by changing appropriate stochastic rate constants. These transitions can be understood in the framework of a simple mathematical model which does not include the detailed biochemical reactions. The model describes a simple gene expression system in which random transitions occur between states 1 and 2, corresponding to low and high levels of protein production. The model has two parameters $`r_1`$ and $`r_2`$ which are the transition rates from state 1 $``$ state 2 and state 2 $``$ state 1 respectively. In the steady state, the probability distribution of protein levels can be calculated. A bimodal distribution corresponds to Type C patterns whereas an unimodal distribution is obtained for Type B patterns. A broad distribution in protein levels (Fig. 14(c)) is obtained corresponding to the gene expression pattern shown in Fig. 12. Transition from one type of distribution to another can be obtained by changing the rate constants $`r_1`$ and $`r_2`$. While experimental evidence for Type A and Type B patterns of gene expression is considerable, we do not know of specific experiments exhibiting Type C patterns. As mentioned in the Introduction, the Type C pattern is similar to that obtained in the case of two-state jump processes in which transitions between two states occur at random time intervals . An example is provided by a spin-$`\frac{1}{2}`$ in the presence of a magnetic field and in contact with a heat bath . The literature on two state jump processes is large and one measurable quantity of interest is the mean first-passage time (MFPT). The time required to switch between two states is a random variable and is known as the first first passage time. Determination of the MFPT and other characteristic measures of the two-state jump processes describing Type C gene expression has not been attempted in this paper. The Type C pattern is reminiscent of a binary digital pulse with states 1 and 2 corresponding to the “0” (OFF) and “1” (ON) states. Genes with expression pattern of Type C may be combined together to construct binary logical circuits. Recently, synthetic networks of genes displaying features of binary logic circuits have been constructed . Bialek  has studied stability and noise in biochemical switches and has shown that switches with long periods of stability and switchability in milliseconds can be constructed from fewer than a hundred molecules. The conclusion is arrived at by studying a model of the synthesis of a single biochemical species, say, proteins in the Langevin formalism. The result obtained is of considerable interest but needs to be verified in a detailed approach involving intermediate processes. The correlation of the amount of random variation in protein distribution, measured by the Fano factor, with the transcriptional and translational rates can be determined in the stochastic simulation approach and the results compared with those obtained in experiments . A detailed study on the dominant contributions to noise in our models of gene expression is in progress and the results will be reported elsewhere.

ACKNOWLEDGEMENTS

I. B. acknowledges helpful discussions with A. M. Kierzek. R. K. was supported by the Council of Scientific and Industrial Research, India under Sanction No. 9/15 (239) / 2002 - EMR - 1.


text:
# Traces of nonextensivity in particle physics due to fluctuations

## 1 Introduction: connection of fluctuations and nonextensivity

Both the notion of fluctuations and that of nonextensivity are nowdays widely known, albeit mostly in the fields of research only indirectly connected with particle physics. Both turns out to be very fruitful and interesting and this is clearly demonstrated by all other lectures given at this workshop (see especially ).

This lecture will be devoted to the case in which evident nonextensivity of some expressions originate in intrinsic fluctuations in the system under consideration (the origin of which is usually not yet fully understood)111Our encounter with this problem is presented in works .. The best introduction to this problem is provided by the observation that in some cosmic ray data (like depth distribution of starting points of cascades in Pamir lead chamber ) one observes clear deviations from the naively expected exponential distributions of some variables which are evidently better fitted by the power-like formulas:

$$\frac{dN}{dh}=\mathrm{const}\mathrm{exp}\left(\frac{h}{\lambda }\right)\mathrm{const}\left[1(1q)\frac{h}{\lambda }\right]^{\frac{1}{1q}}.$$
(1)
Here $`N`$ denotes the number of counts at depth $`h`$ (cf.  for details). Whereas in  we have proposed as explanation a possible fluctuations of the mean free path $`\lambda `$ in eq. (1) characterised by relative variance $`\omega =\frac{\left(\sigma ^2\sigma ^2\right)}{\sigma ^2}\mathrm{\hspace{0.17em}0.2}`$, in  the same data were fitted by power-like (Lévy type) formula as above keeping $`\lambda `$ fixed and setting $`q=1.3`$. In this way we have learned about Tsallis statistics and Tsallis nonextensive entropy and distributions222See Tsallis lecture  and references therein (cf. also ) for detailed necessary information concerning Tsallis statistics and the non-extensivity.. By closer inspection of the above example we have been able to propose a new physical meaning of the nonextensivity parameter $`q`$, as a measure of intrinsic fluctuations in the system . Fluctuations are therefore proposed as a new source of nonextensivity which should be added  to the previously known and listed in literature sources of the nonextensivity (like long-range correlations, memory effects or fractal-like structure of the corresponding phases space ) .

To demonstrate this conjecture let us notice that for $`q>1`$ case, where $`\epsilon (0,\mathrm{})`$, one can write a kind of Mellin transform (here $`\alpha =\frac{1}{q1}`$) :

$`\left(1+{\displaystyle \frac{\epsilon }{\lambda _0}}{\displaystyle \frac{1}{\alpha }}\right)^a`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Gamma }(\alpha )}}{\displaystyle _0^{\mathrm{}}}𝑑\eta \eta ^{\alpha 1}\mathrm{exp}\left[\eta \left(1+{\displaystyle \frac{\epsilon }{\lambda _0}}{\displaystyle \frac{1}{\alpha }}\right)\right]`$ (2)
$`=`$ $`{\displaystyle _0^{\mathrm{}}}\mathrm{exp}\left({\displaystyle \frac{\epsilon }{\lambda }}\right)f_{q>1}\left({\displaystyle \frac{1}{\lambda }}\right)d\left({\displaystyle \frac{1}{\lambda }}\right),`$
where $`f_{q>1}(1/\lambda )`$ is given by the following gamma distribution:

$$f_{q>1}\left(\frac{1}{\lambda }\right)=f_{q>1}(\frac{1}{\lambda },\frac{1}{\lambda _0})=\frac{\mu }{\mathrm{\Gamma }(\alpha )}\left(\frac{\mu }{\lambda }\right)^{\alpha 1}\mathrm{exp}\left(\frac{\mu }{\lambda }\right)$$
(3)
with $`\mu =\alpha \lambda _0`$ and with mean value and variation in the form:

$$\frac{1}{\lambda }=\frac{1}{\lambda _0}\mathrm{and}\left(\frac{1}{\lambda }\right)^2\frac{1}{\lambda }^2=\frac{1}{\alpha \lambda _0^2}.$$
(4)
For the $`q<1`$ case $`\epsilon `$ is limited to $`\epsilon [0,\lambda _0/(1q)]`$. Proceeding in the same way as before (but with $`\alpha ^{}=\alpha =\frac{1}{1q}`$) one gets:

$`\left[1{\displaystyle \frac{\epsilon }{\alpha ^{}\lambda _0}}\right]^\alpha ^{}`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Gamma }(\alpha ^{})}}{\displaystyle _0^{\mathrm{}}}𝑑\eta \eta ^{\alpha ^{}1}\mathrm{exp}\left[\eta \left(1+{\displaystyle \frac{\epsilon }{\alpha ^{}\lambda _0\epsilon }}\right)\right]`$ (5)
$`=`$ $`{\displaystyle _0^{\mathrm{}}}\mathrm{exp}\left({\displaystyle \frac{\epsilon }{\lambda }}\right)f_{q<1}\left({\displaystyle \frac{1}{\lambda }}\right)d\left({\displaystyle \frac{1}{\lambda }}\right),`$
where $`f_{q<1}(1/\lambda )`$ is given by the same gamma distribution as in (3) but this time with $`\alpha \alpha ^{}`$ and $`\mu =\mu (\epsilon )=\alpha ^{}\lambda _0\epsilon `$. Contrary to the $`q>1`$ case, this time the fluctuations depend on the value of the variable in question, i.e., the mean value and variance are now both $`\epsilon `$-dependent:

$$\frac{1}{\lambda }=\frac{1}{\lambda _0\frac{\epsilon }{\alpha ^{}}}\mathrm{and}\left(\frac{1}{\lambda }\right)^2\frac{1}{\lambda }^2=\frac{1}{\alpha ^{}}\frac{1}{\left(\lambda _0\frac{\epsilon }{\alpha ^{}}\right)^2}.$$
(6)
However, in both cases the relative variances,

$$\omega =\frac{\left(\frac{1}{\lambda }\right)^2\frac{1}{\lambda }^2}{\frac{1}{\lambda }^2}=\left\{\begin{array}{c}\frac{1}{\alpha }\\ \\ \frac{1}{\alpha ^{}}\end{array}\right\}=\left\{\begin{array}{c}q1\\ \\ 1q\end{array}\right\},$$
(7)
remain $`\epsilon `$-independent and depend only on the parameter $`q`$ 333Notice that, with increasing $`\alpha `$ or $`\alpha ^{}`$ (i.e., for $`q1`$) both variances (7) decrease and asymptotically gamma distributions (3) becomes a delta function, $`f(1/\lambda )=\delta (\lambda \lambda _0)`$.. It means therefore that  (at least for the fluctuations distributed according to gamma distribution)

$$L=\mathrm{exp}\left(\frac{x}{\lambda _0}\right)L_q=\mathrm{exp}_q\left(\frac{x}{\lambda _0}\right)=\mathrm{exp}\left(\frac{x}{\lambda }\right),$$
(8)
with $`q=1\pm \omega `$ for $`q>1(+)`$ and $`q<1()`$, i.e., there is connection between the measure of fluctuations $`\omega `$ and the measure of nonextensivity $`q`$ (it has been confirmed recently in ).

## 2 Where are the fluctuations coming from?

### 2.1 Generalities

The proposed interpretation of the parameter $`q`$ leads immediately to the following question: why and under what circumstances it is the gamma distribution that describes fluctuations of the parameter $`\lambda `$? To address this question let us write the usual Langevin equation for the stochastic variable $`\lambda `$ :

$$\frac{d\lambda }{dt}+\left[\frac{1}{\tau }+\xi (t)\right]\lambda =\varphi =\mathrm{const}>\mathrm{\hspace{0.17em}0}.$$
(9)
with damping constant $`\tau `$ and with source term $`\varphi `$, different for the two cases considered, namely:

$$\varphi =\varphi _{q<1}=\frac{1}{\tau }\left(\lambda _0\frac{\epsilon }{\alpha ^{}}\right)\mathrm{whereas}\varphi =\varphi _{q>1}=\frac{\lambda _0}{\tau }.$$
(10)
For the usual stochastic processes defined by the white gaussian noise form of $`\xi (t)`$ 444It means that ensemble mean $`\xi (t)=\mathrm{\hspace{0.17em}0}`$ and correlator (for sufficiently fast changes) $`\xi (t)\xi (t+\mathrm{\Delta }t)=\mathrm{\hspace{0.17em}2}D\delta (\mathrm{\Delta }t)`$. Constants $`\tau `$ and $`D`$ define, respectively, the mean time for changes and their variance by means of the following conditions: $`\lambda (t)=\lambda _0\mathrm{exp}\left(\frac{t}{\tau }\right)`$ and $`\lambda ^2(t=\mathrm{})=\frac{1}{2}D\tau `$. Thermodynamical equilibrium is assumed here (i.e., $`t>>\tau `$, in which case the influence of the initial condition vanishes and the mean squared of $`\lambda `$ has value corresponding to the state of equilibrium). one obtains the following Fokker-Plank equation for the distribution function of the variable $`\lambda `$

$$\frac{df(\lambda )}{dt}=\frac{}{\lambda }K_1f(\lambda )+\frac{1}{2}\frac{^2}{\lambda ^2}K_2f(\lambda ),$$
(11)
where the intensity coefficients $`K_{1,2}`$ are defined by eq.(9) and are equal to :

$$K_1(\lambda )=\varphi \frac{\lambda }{\tau }+D\lambda \mathrm{and}K_2(\lambda )=\mathrm{\hspace{0.17em}2}D\lambda ^2.$$
(12)
From it we get the following expression for the distribution function of the variable $`\lambda `$:

$$f(\lambda )=\frac{c}{K_2(\lambda )}\mathrm{exp}\left[\mathrm{\hspace{0.17em}2}_0^\lambda 𝑑\lambda ^{}\frac{K_1(\lambda ^{})}{K_2(\lambda ^{})}\right]$$
(13)
which is, indeed, a gamma distribution (3) in variable $`1/\lambda `$, with the constant $`c`$ defined by the normalization condition, $`_0^{\mathrm{}}d(1/\lambda )f(1/\lambda )=1`$, and depending on two parameters: $`\mu (\epsilon )=\frac{\varphi _q(\epsilon )}{D}`$ and $`\alpha _q=\frac{1}{\tau D}`$ with $`\varphi _q=\varphi _{q>1,q<1}`$ and $`\alpha _q=(\alpha ,\alpha ^{})`$ for, respectively, $`q>1`$ and $`q<1`$. This means that we have obtained eqs.(7) with $`\omega =\frac{1}{\tau D}`$ and, therefore, the parameter of nonextensivity $`q`$ is given by the parameter $`D`$ and by the damping constant $`\tau `$ describing the white noise.

### 2.2 Temperature fluctuations

The above discussion rests on the stochastic equation (9). Therefore the previously asked question is not yet fully answered but can be refrazed in the following way: can one point the possible physical situation where such fluctuations could be present in the realm of particle physics? Our proposition to answer it is to identify $`\lambda =T`$, i.e., to concentrate on the possibility of fluctuations of temperature widely discussed in the literature in different contexts . In all cases of interest to us temperature $`T`$ is variable encountered in statistical descriptions of multiparticle collision processes . Our reasoning goes then as follows: Suppose that we have a thermodynamic system, in a small (mentally separated) part of which the temperature fluctuates with $`\mathrm{\Delta }TT`$. Let $`\lambda (t)`$ describe stochastic changes of the temperature in time. If the mean temperature of the system is $`T=T_0`$ then, as result of fluctuations in some small selected region, the actual temperature equals $`T^{}=T_0\tau \xi (t)T`$ and the inevitable exchange of heat between this selected region and the rest of the system leads to the process of equilibration of the temperature which is described by the following equation of the type of eq. (9)  :

$$\frac{T}{t}\frac{1}{\tau }(T^{}T)+\mathrm{\Omega }_q=0$$
(14)
(here $`\mathrm{\Omega }_{q<1}=\frac{\epsilon }{\tau \alpha ^{}}`$ and $`\mathrm{\Omega }_{q>1}=0`$). In this way we have demonstrated that eq. (9) can, indeed, be used in statistical models of multiparticle production as far as they allow the $`q`$-statistic as their basis.

The $`q<1`$ case requires a few additional words of explanation. Notice the presence of the internal heat source, which is dissipating (or transfering) energy from the region where (due to fluctuations) the temperature $`T`$ is higher to colder ones. It could be any kind of convection-type flow of energy; for example, it could be connected with emission of particles from that region. The heat release given by $`\epsilon /(\tau \alpha ^{})`$ depends on $`\epsilon `$ (but it is only a part of $`\epsilon `$ that is released). In the case of such energy release (connected with emission of particles) there is additional cooling of the whole system. If this process is sufficiently fast, it could happen that there is no way to reach a stationary distribution of temperature (because the transfer of heat from the outside can be not sufficient for the development of the state of equilibrium). On the other hand (albeit this is not our case here) in the reverse process we could face the ”heat explosion” situation (which could happen if the velocity of the exothermic burning reaction grows sufficiently fast; in this case because of nonexistence of stationary distribution we have fast nonstationary heating of the substance and acceleration of the respective reaction)555It should be noticed that in the case of $`q<1`$ temperature does not reach stationary state because, cf. Eq. (6), $`1/T=\mathrm{\hspace{0.17em}1}/(T_0\epsilon /\alpha ^{})`$, whereas for $`q>1`$ we had $`<1/T>=1/T_0`$. As a consequence the corresponding Lévy distributions are defined only for $`\epsilon (0,T_0\alpha ^{}`$) because for $`\epsilon T_0\alpha ^{}`$, $`<T>0`$. Such asymptotic (i.e., for $`t/\tau \mathrm{}`$) cooling of the system ($`T0`$) can be also deduced form Eq. (14) for $`\epsilon T_0\alpha ^{}`$..

The most interesting example of the possible existence of such fluctuations is the trace of power like behaviour of the transverse momentum distribution in multiparticle production processes encountered in heavy ion collisions . Such collisions are of special interest because they are the only place where new state of matter, the Quark Gluon Plasma, can be produced . Transverse momentum distributions are believed to provide information on the temperature $`T`$ of reaction, which is given by the inverse slope of $`dN/dp_T`$, if it is exponential one. If it is not, question arises what we are really measuring. One explanation is the possible flow of the matter, the other, which we shall follow here, is the nonextensivity (or rather fluctuations leading to it). Namely, as was discussed in detail in  the extreme conditions of high density and temperature occuring in ultrarelativistic heavy ion collisions can invalidate the usual BG approach and lead to $`q>1`$, i.e., to

$$\frac{dN(p_T)}{dp_T}=\mathrm{const}\left[1(1q)\frac{\sqrt{m^2+p_T^2}}{T}\right]^{\frac{1}{1q}}.$$
(15)
Here $`m`$ is the mass of produced particle and $`T`$ is, for the $`q=1`$ case, the temperature of the hadronic system produced. Although very small ($`|q1|0.015`$) this deviation, if interpreted according to eq. (8)), leads to quite large relative fluctuations of temperature existing in the nuclear collisions, $`\mathrm{\Delta }T/T0.12`$. It is important to stress that these are fluctuations existing in small parts of hadronic system in respect to the whole system rather than of the event-by-event type for which, $`\mathrm{\Delta }T/T=0.06/\sqrt{N}0`$ for large $`N`$ (cf.  for relevant references). Such fluctuations are potentially very interesting because they provide direct measure of the total heat capacity $`C`$ of the system,

$$\frac{\sigma ^2(\beta )}{\beta ^2}=\frac{1}{C}=\omega =q1,$$
(16)
($`\beta =\frac{1}{T}`$) in terms of $`\omega =q1`$. Therefore, measuring both the temperature of reaction $`T`$ and (via nonextensivity $`q1`$) its total heat capacity $`C`$, one can not only check whether an approximate thermodynamics state is formed in a single collision but also what are its theromdynamical properties (especially in what concerns the existence and type of the possible phase transitions ).

To observe such fluctuations an event-by-event analysis of data is needed . Two scenarios must be checked: $`(a)`$ $`T`$ is constant in each event but because of different initial conditions it fluctuates from event to event and $`(b)`$ $`T`$ fluctuates in each event around some mean value $`T_0`$. Fig. 1 shows typical event obtained in simulations performed for central $`Pb+Pb`$ collisions taking place for beam energy equal $`E_{beam}=3A`$TeV in which density of particles in central region (defined by rapidity window $`1.5<y<1.5`$) is equal to $`\frac{dN}{dy}=6000`$ (this is the usual value given by commonly used event generators ). In case $`(a)`$ in each event one expects exponential dependence with $`T=T_{event}`$ and possible departure from it would occur only after averaging over all events. It would reflect fluctuations originating from different initial conditions for each particular collision. This situation is illustrated in Fig. 1a where $`p_T`$ distributions for $`T=200`$ MeV (black symbols) and $`T=250`$ MeV (open symbols) are presented. Such values of $`T`$ correspond to typical uncertainties in $`T`$ expected at LHC accelerator at CERN. Notice that both curves presented here are straight lines. In case $`(b)`$ one should observe departure from the exponential behaviour already on the single event level and it should be fully given by $`q>1`$. It reflects situation when, due to some intrinsically dynamical reasons, different parts of a given event can have different temperatures . In Fig. 1b black symbols represent exponential dependence obtained for $`T=200`$ MeV (the same as in Fig. 1a), open symbols show the power-like dependence as given by (15) with the same $`T`$ and with $`q=1.05`$ (notice that the corresponding curve bends slightly upward here). In this typical event we have $`18000`$ secondaries, i.e., practically the maximal possible number. Notice that points with highest $`p_T`$ correspond already to single particles. As one can see, experimental differentiation between these two scenarios will be very difficult, although not totally impossible. On the other hand, if successful it would be very rewarding - as we have stressed before.

One should mention at this point that to the same cathegory of fluctuating temperature belongs also attempt  to fit energy spectra in both the longitudinal and transverse momenta of particles produced in the $`e^+e^{}`$ annihilation processes at high energies, novel nonextensive formulation of Hagedorn statistical model of hadronization process  and description of single particle spectra .

### 2.3 Nonexponential decays

Another hint for intrinsic fluctuations operating in the physical system could be the known phenomenon of nonexponential decays . Spontaneous decays of quantum-mechanical unstable systems cannot be described by the pure exponential law (neither for short nor for long times) and survival time probability is $`P(t)t^\delta `$ instead of exponential one. It turns out  that by using random matrix approach, such decays can emerge in a natural way from the possible fluctuations of parameter $`\gamma =1/\tau `$ in the exponential distribution $`P(t)=\mathrm{exp}(\gamma t)`$. Namely, in the case of multichannel decays (with $`\nu `$ channels of equal widths involved) one gets fluctuating widths distributed according to gamma function

$$P_\nu (\gamma )=\frac{1}{\mathrm{\Gamma }\left(\frac{\nu }{2}\right)}\left(\frac{\nu }{2<\gamma >}\right)\left(\frac{\nu \gamma }{2<\gamma >}\right)^{\frac{\nu }{2}1}\mathrm{exp}\left(\frac{\nu \gamma }{2<\gamma >}\right)$$
(17)
and strength of their fluctuations is given by relative variance $`\frac{(\gamma <\gamma >)^2}{<\gamma >^2}=\frac{2}{\nu }`$, which decreases with increasing $`\nu `$. According to , it means therefore that, indeed,

$$L_q(t,\tau _0)=\frac{2q}{\tau _0}\left[1(1q)\frac{t}{\tau _0}\right]^{\frac{1}{1q}},$$
(18)
with the nonextensivity parameter equal to $`q=1+\frac{2}{\nu }`$.

## 3 Summary

There is steadily growing evidence that some peculiar features observed in particle and nuclear physics (including cosmic rays) can be most consistently explained in terms of the suitable applications of nonextensive statistic of Tsallis. Here we were able to show only some selected examples, more can be found in . However, there is also some resistance towards this idea, the best example of which is provided in . It is shown there that mean multiplicity of neutral mesons produced in $`p\overline{p}`$ collisions as a function of their mass (in the range from $`m_\eta =0.55`$ GeV to $`M_\mathrm{{\rm Y}}=9.5`$ GeV) and the transverse mass $`m_T`$ spectra of pions (in the range of $`m_T1÷15`$ GeV), both show a remarkable universal behaviour following over $`10`$ orders of magnitude the same power law function $`Cx^P`$ (with $`x=m`$ or $`x=m_T`$) with $`P10.1`$ and $`P9.6`$, respectively. In this work such a form was just postulated whereas it emerges naturally in $`q`$-statistics with $`q=1+1/P1.1`$ (quite close to results of ). We regard it as new, strong imprint of nonextensivity present in multiparticle production processes (the origin of which remains, however, to be yet discovered). This interpretation is additionally supported by the fact that in both cases considered in  the constant $`c`$ is the same. Apparently there is no such phenomenon in $`AA`$ collisions which has simple and interesting explanation: in nuclear collisions volume of interaction is much bigger what makes the heat capacity $`C`$ also bigger. This in turn, cf. eq.(16), makes $`q`$ smaller. On should then, indeed, expect that $`q_{hadronic}>>q_{nuclear}`$, as observed.

As closing remark let us point out the alternate way of getting nonextensive (i.e., with $`q1`$) distributions for thermal models (cf. our remarks in  and the more recent ideas presented in ). Notice that if we allow for the temperature $`T`$ to be energy $`E`$ dependent, i.e., that $`T=T(E)=T_0+a(E\stackrel{~}{E})`$ (with $`a=1/C_V`$) then the usual equation on the probability $`P(E)`$ that a system $`A`$ (interacting with the heat bath $`A^{}`$ with temperature $`T`$) has energy $`E`$,

$$d\mathrm{ln}[P(E)]dE/TP(E)\mathrm{exp}\left(\frac{E}{T}\right)$$
(19)
becomes

$$d\mathrm{ln}P(E)\frac{1}{T_0+a(E\stackrel{~}{E})}dEP(E)\left[1(1q)\frac{E\stackrel{~}{E}}{T_0}\right]^{\frac{1}{1q}}$$
(20)
with $`q=1+a`$. This approach could then find its possible application to studies of fluctuations on event-by-event basis  (with all reservations expressed in  accounted for).

## Acknowledgments

GW would like to thank Prof. N.G. Antoniou and all Organizers of X-th International Workshop on Multiparticle Production, Correlations and Fluctuations in QCD for financial support and kind hospitality.