text:
# Biological Networks

## 1 Introduction

Networks are widely prevalent in all spheres of life . A network of acquaintances is the simplest example one can think of. Social, economic and political networks of various kinds are part of human society. The internet, a network of information resources, plays a vital role in the gathering, sharing and transmission of information. A network consists of nodes connected by links. Figure 1 shows the example of a network in which the solid circles denote the nodes and the solid lines the links. Some examples of real life networks are as follows: in a network describing an electrical power grid, the generators, transformers and substations are the nodes and the high-voltage transmission lines connecting them the links. In the World Wide Web (WWW), the documents/pages constitute the nodes. These are connected to other documents/pages through links. In a collaboration graph of movie actors, the nodes represent the actors. Two actors are connected by a link if they appear in the same movie. In a citation network, the nodes are the papers published in refereed journals. A paper is linked to all the other papers it cites. Cellular processes are controlled by various types of biochemical networks. A metabolic network  controls the processes which generate mass and energy from nutritive matter. The nodes in such a network are the substrates such as ATP, ADP and $`H_2O`$. Two substrates are connected by a link if both of them participate in the same biochemical reaction. Traditional cell biology assigns specific functional roles to individual proteins, such as catalysts, signalling molecules and constituents of cellular matter. In the post-genomic era, there is an increasing emphasis on understanding the functions of proteins as parts of an interacting network and also on the collective, emergent properties of the network. In a protein-protein interaction network , the nodes represent the proteins. A link exists between two nodes if the corresponding proteins have a direct physical interaction.

The networks discussed above have complex topology. Spectacular advances in computerisation of data acquisition (the Human Genome Project is a prime example) have made it possible to construct large databases which contain information on the topology of real life networks. The advent of powerful computers has given rise to extensive investigations of networks containing millions of nodes. The interesting fact emerging out of these studies is that biological networks share common topological features with non-biological networks. There appears to be a general blueprint for the large scale organisation of several of these networks. In Section 2 of this review, we discuss two types of biological networks, namely, the metabolic networks of several organisms and the protein-protein interaction networks associated with the yeast S. cerevisiae  and the human gastric pathogen H. Pylori . The major topological features of these networks are described and the similarity in the design principles of large-scale biological and non-biological networks pointed out.

Gene regulatory networks are the most significant examples of biological networks. Gene expression and regulation are the central activities of a living cell . Genes are fragments of $`DNA`$ molecules and determine the structure of functional molecules like $`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. The information present in the gene is expressed through the processes of transcription and translation. During transcription, the sequence along one of the strands of the $`DNA`$ molecule is copied onto a $`RNA`$ molecule ($`mRNA`$ ). The sequence of the $`mRNA`$ molecule is then translated into the sequence of amino acids, which determines the functional nature of the protein molecule produced. In a gene regulatory network, the protein encoded by one gene can regulate the expression of other genes. These genes in turn produce new regulatory proteins which control still other genes. A protein may also regulate its own level of production through an autoregulatory feedback process. The occurrence of cell differentiation, when an organism grows from its embryonic stage, depends upon the selective switching on of gene expression in individual cells. All these cells have identical sets of genes but follow different developmental pathways depending upon the patterns of gene expression in the cells. Thus distinct types of cells such as hair and skin cells are obtained. Gene expression is also regulated in metabolism and progression through the cell cycle as well as in responses to external signals. Infected cells can multiply because the expression of certain genes is “on” in these cells whereas in normal cells the expression of the same genes is “off”.

Despite a vast amount of experimental data, the complex dynamical processes involved in gene regulation are not fully understood as yet. A large number of theoretical studies has been undertaken  but only a few of these make quantitative predictions in agreement with experimental results. Two key concepts which emerge out of the theoretical studies are: nonlinearity of the network dynamics and the role of stochasticity in gene expression and regulation . The variables of interest in the network dynamics are the concentrations of the $`mRNAs`$, proteins and other biomolecules within the cell. The rate of change in the concentration of a biomolecular species is a nonlinear function of the other variables. The dynamics is governed by a set of coupled non-linear differential equations which in most cases are solved numerically. Let $`U`$ be the concentration of, say, a particular type of protein in the cell. The rate of change of $`U`$ is given by

$`\frac{dU}{dt}=`$ (Production - Loss/Decay ) of $`U`$ per unit time.
The production term is a nonlinear function of the other concentration variables and the loss term is usually proportional to $`U`$. In Section 3 of this review, we briefly describe the major features of nonlinearity in the dynamics of gene networks. There is currently a significant emerging trend to utilise the concepts and techniques of nonlinear physics in the actual construction of synthetic gene regulatory networks with a variety of applications. In Section 3, a specific example of this, namely, the genetic toggle switch  will be given.

The biochemical rate equations which govern the dynamics of gene regulatory networks are deterministic in nature. Many molecules associated with the networks have low intracellular concentrations and consequently fluctuations in reaction rates are considerably large. Gene expression involves a series of biochemical reactions and due to stochastic fluctuations in the reaction rates, proteins are produced in short bursts at random time intervals. In the last few years, there is an increasing realization that stochasticity plays a significant role in biological processes . To give an example, consider the situation in which two independently produced regulatory proteins A and B are in competition to control a developmental switch that selects between two pathways depending on which protein wins. The protein concentrations have to reach effective levels in order to activate the switch. Due to stochastic fluctuations, the amounts of proteins A and B produced as a function of time can vary widely from cell to cell. In some cells, protein A reaches the effective level first and activates the developmental switch along one specific pathway. In the other cells, protein B takes control and the other pathway is activated. Thus even a clonal cell population exhibits phenotypic variations as the cells follow different developmental pathways. Environmental signals can bias the probabilities of path choice in a regulatory circuit. Organisms make use of this mechanism to increase the probability of survival in a hostile environment. Cells often utilise fluctuations (noise) to randomize choices of developmental pathways when such randomization is desirable for the survival and growth of the organism. A well-known example is that of the phage $`\lambda `$ lysis-lysogeny network . The bacterial E.coli cells, when infected by the virus phage $`\lambda `$ , can follow two developmental pathways: lysis and lysogeny. In the lysogenic state, the infection is dormant. Phage $`\lambda `$ is inert and integrated into the host cell’s chromosome. It replicates along with the bacterial DNA and each new cell contains the dormant phage. In the lytic state, the infection proliferates. The viral DNA replicates using the host cell machinery giving rise to a large number of progeny phage. These in turn lyse or burst the host bacterium cell and the infection spreads to more cells. Again, due to stochastic fluctuations, the cell population divides into two subpopulations: lysogenic and lytic. The selection of a developmental pathway after the host cell is infected is not deterministic but probabilistic. In Section 4 of this review, the effect of stochasticity on the dynamics of the $`\lambda `$-phage network will be briefly discussed. The network illustrates the competitive control of a developmental switch by two regulatory proteins.

As already mentioned before, gene expression/regulation involves several biochemical reactions with appreciable stochastic fluctuations. Gillespie  has proposed a Monte Carlo simulation algorithm to describe the kinetics of coupled stochastic reactions. This method is physically more rigorous than the conventional differential equation approach. The inherent assumption in the latter method is that the temporal changes in the concentrations of reacting molecules are both continuous and deterministic. The assumption is not true if the concentrations are small and the reaction rates slow or if the system undergoes large, rapid and discrete transitions. In the Gillespie algorithm, changes in the numbers of the reacting molecules occur in integral numbers brought about by random, distinct reaction events. The Gillespie algorithm is described in detail in Section 4 and some illustrative examples are given. Recent experiments at the level of a single cell have shown that gene expression occurs in abrupt stochastic bursts . Further, in an ensemble of cells, the levels of proteins produced have a bimodal distribution. In a large fraction of cells, the gene expression is either off or has a high value. We have proposed a model of gene expression the essential features of which are stochasticity and cooperative binding of RNA polymerase, the molecule responsible for transcription . The model can reproduce the bimodality observed in experiments. We include a description of the model in Section 4 to give an additional example of the effect of stochasticity on gene regulation. Section 5 of the review contains concluding remarks.

The emphasis in this review is on recent studies of biological networks. There are a number of exhaustive reviews and books on earlier work. The three major themes that several recent studies focus on are: topology of networks, nonlinear dynamics and its consequences and the role of stochasticity in biological processes. The present review is meant to be an introduction to these themes and to highlight the fact that interdisciplinary approaches are essential to develop an integrated understanding of biological networks.

## 2 Topology of complex networks

Many real life networks have a complex structure. The mathematicians Erdös and Rényi  were the first to propose a model of a complex network known as a random graph. One starts with $`N`$ nodes and connects every pair of nodes with probability $`p`$. The graph thus has approximately $`\frac{pN(N1)}{2}`$ links distributed in a random manner. Studies of real life networks, however, reveal that these cannot be described as random graphs. This distinction is possible on the basis of quantitative measurements of certain topological features which we define below. Several complex networks including the random graph are described as small world networks . The small world idea implies that though the networks are large in size (the number of nodes in a network is a measure of its size), any pair of nodes can be connected by a short path. The distance between two nodes is given by the number of links along the shortest path connecting the nodes. In Figure 1, the distance between the nodes A and B is three. The diameter of the network, also known as the average path length l, is the average of the distances between all pairs of nodes. The global population is huge but still we live in a small world as any random pair of individuals are connected to each other through a short path of intermediate acquaintances. This was first established by Stanley Milgram  who found out that the average path length of intermediate acquaintances is six. In a small world network, the diameter scales as the logarithm of the number of nodes.

The second measurable topological feature of a complex network is its degree distribution . The number of links by which a node is connected to the other nodes varies from node to node. Let $`P(k)`$ be the probability that a randomly selected node has exactly $`k`$ links. Equivalently, $`P(k)`$ is the fraction of nodes, on an average, which has exactly $`k`$ links. One can define an average degree $`k`$ of the network, the degree of a node being the number of links attached to the node. In a random graph, the links are established randomly and most of the nodes have degrees close to $`k`$ . The degree distribution $`P(k)`$ vs. $`k`$ is Poissonian. It is strongly peaked at $`k=k`$ and has an exponential decay for large k, i.e., $`P(k)e^k`$ for $`kk`$ and $`kk`$ . In many real life networks, the degree distribution $`P(k)`$ has no well-defined peak but has a power-law distribution

$$P(k)k^\gamma $$
(1)
where the exponent $`\gamma `$ is a numerical constant. Such networks are known as scale-free networks because they are not tied to a specific scale. $`P(k)`$ has a finite value over a wide range of $`k`$ values. The power-law form of the degree distribution implies that the networks are extremely inhomogeneous unlike in the case of a random graph. In a scale- free network, there are many nodes with few links and a few nodes with many links. The highly connected nodes play a key role in the functionality of the network. Both the random graph and the scale-free networks are small world networks.

The third topological quantity which is measurable is known as the clustering coefficient . The coefficient is a measure of the tendency of the nodes of the network towards clustering. In a social network, the individuals are the nodes and two nodes are connected by a link if the individuals are acquainted with each other. In such a network, one’s friend’s friends are also likely to be one’s friends giving rise to a clustering of acquaintances. The clustering coefficient is defined in the following manner. Let us select a specific node $`i`$ in the network which is connected by $`k_i`$ links to $`k_i`$ other nodes. If these first neighbours are all part of a cluster, there would be $`\frac{k_i(k_i1)}{2}`$ links between them. The clustering coefficient $`C_i`$ of node i is given by

$$C_i=\frac{2E_i}{k_i(k_i1)}$$
(2)
where $`E_i`$ is the number of actual links which exist between the $`k_i`$ nodes. The clustering coefficient C of the whole network is obtained by taking an average over all the $`C_i`$ values. The utility of the clustering coefficient is demonstrated in the following example. The neural network of the nematode worm C. Elegans is small in size . The number of neurons which constitute the nodes of the network is 282. A link exists between two nodes if the neurons are connected by either a synapse or a gap junction. The average degree of the network is $`k=`$14. Now consider a random graph of the same size and average degree. The average path lengths for the neural network and the random graph are similar, 2.65 and 2.25 respectively. Is the neural network then a random graph? The answer is no as the clustering coefficient of the former has the value 0.28 which is much larger than the value 0.05 in the case of the latter network. Examples of real life networks which are scale-free are : the collaboration graph of movie actors (size $`N`$ of the network = 212 250 nodes, average degree $`k=`$ 28.78, the exponent $`\gamma `$ in Eq.(1) is $`\gamma `$ = 2.3 ), the WWW ($`N=`$ 325 729, $`k=`$ 5.46, $`\gamma =`$ 2.1) and the network of citations ($`N=`$ 783 339 papers, $`k=`$ 8.57, $`\gamma =`$ 3). The results are obtained from available databases. A more comprehensive and up to date list of networks is given in Ref. .

We now discuss some complex, biological networks. Recently, Jeong et al  have systematically investigated the topological properties of the core metabolic networks of 43 different organisms representing all the three domains of life. The data on these organisms are available in the WIT (What Is There) database. As already mentioned in the Introduction, the nodes of the metabolic network are the different substrates. Two substrates are connected by a link if they participate in the same biochemical reaction. The metabolic networks have different sizes, the less complex organisms having smaller sizes. There is considerable variation in the individual constituents and the pathways of the networks. Yet they display identical topological scaling properties which resemble those of complex non-biological networks. The metabolic networks have been found to belong to the class of scale-free networks. The probability that a substrate participates in k reactions has a power-law distribution. The links in a metabolic network are directed as many biochemical reactions are preferentially catalysed in one direction. For each node, one has to distinguish between incoming and outgoing links. Correspondingly, there are two exponents $`\gamma _{in}`$ and $`\gamma _{out}`$ . The exponents turn out to have the same value of 2.2.

In the metabolic network, the distance between two substrates is given by the number of links (reactions) in the shortest biochemical pathway connecting the two substrates. A surprising result obtained by Jeong et al is that the diameter of the metabolic network is the same for all the 43 organisms, i.e., it does not depend upon the number of substrates (nodes) belonging to the network. This is counterintuitive and only possible if with increasing organism complexity pre-existing individual substrates are increasingly connected in order to maintain a more or less constant network diameter. In support of this conjecture, Jeong et al found that the average number of reactions in which a certain substrate participates increases as the number of substrates in the organism increases. Conservation of the network diameter may be favourable for the survival and growth of an organism. A larger diameter would possibly diminish the organism’s ability to respond to changes in an efficient manner.

The scale-free character of the metabolic network implies that a few hubs which are highly connected play a dominant role in the functioning of the network. On sequential removal of these highly-connected nodes, the network diameter rises sharply and ultimately the network disintegrates into isolated fragments. On the other hand, the network diameter does not change appreciably when the nodes with a few links are removed from the network. Scale-free networks, in general, are robust against random mutations/errors but vulnerable to attacks targeted at highly connected nodes. Complex communication networks are surprisingly robust, local failures rarely hamper the global transmission of information. Organisms can grow and survive in hostile environments due to the error tolerance of the underlying metabolic network. A random graph is not as robust against random mutations/errors. Mutagenesis studies in-silico and in-vivo  have established the remarkable error tolerance of the metabolic network of E.coli on removing a large number of metabolic enzymes. Jeong et al, in their study of the metabolic networks of organisms found that only $`4\%`$ of all the substrates present in all the 43 organisms are present in all the species. The striking fact is that the small number of substrates, common to all species, turn out to be the most highly connected ones. On the other hand, there are species-specific differences in the case of less-connected substrates.

Jeong et al in a separate study  have investigated the protein-protein interaction network of the yeast S.cerevisiae. The network has 1870 proteins as nodes which are linked by 2240 direct physical interactions identified mostly by systematic two-hybrid experiments. Actual measurements show that the probability $`P(k)`$ that a given yeast protein interacts with $`k`$ other yeast proteins has a power-law distribution with an exponential cutoff at $`k_c=`$ 20.

$$P(k)(k+k_0)^\gamma e^{\frac{(k+k_0)}{k_c}}$$
(3)
with $`k_0=`$ 1 and $`\gamma =`$ 2.4. The protein-protein interaction network of the bacterium H. Pylori  displays similar topology. For the metabolic networks, the exponent $`\gamma `$ has the value 2.2. The value of $`\gamma `$ falls in the range 2.0-2.5 for many scale-free networks. Like the metabolic network and other scale-free networks, the protein-protein interaction network is found to be immune to random mutations. The removal of highly connected nodes may, however, disrupt the network function. The protein product of the p53 tumor-suppressor gene is one of the most highly connected proteins found in human cells. Mutations of p53 gene affect cellular functions severely from a biomedical point of view.

In fact, Jeong et al’s study on the protein-protein interaction network in yeast shows that proteins with five or lesser number of links constitute $``$ 93% of the total number of proteins but only $``$ 21% of them are essential so that the removal of such proteins proves to be lethal. In contrast, only $``$ 0.7% of the total number of proteins have more than 15 links but single deletion of $``$ 62% of these severely affects the functioning of the network. It is possible that the proteins which constitute the highly connected nodes in a network share common structural features. These features favour the binding of many different types of proteins to the proteins in question. The scale-free character of both the metabolic and protein-protein interaction networks suggests the evolutionary selection of a common large scale structure of biological networks. Studies of other biological networks are expected to provide further evidence for this idea.

## 3 Nonlinear dynamics

The dynamics of gene regulatory networks are described by coupled nonlinear p.d.e.’s which can be collectively represented as

$$\frac{dX(t)}{dt}=f(X,R)$$
(4)
where $`X(t)`$ is the $`N`$ -component state vector $`(X_1(t),\mathrm{},X_N(t))`$ and $`f`$ is a set of nonlinear functions $`f_1(X,R),\mathrm{}.,f_N(X,R)`$ . There are thus N coupled p.d.e.’s and an individual p.d.e. is of the form

$$\frac{dX_i(t)}{dt}=f_i(X_1,\mathrm{}.X_N,R)$$
(5)
There are in total N species of biochemical molecules participating in M reactions. $`X_i(t)`$ $`(i=1,\mathrm{},N)`$ represents the concentration of the ith molecular species at time $`t`$ . R represents a set of control parameters. The functions $`f_i^{}s`$ are nonlinear functions of the $`X_i^{}s`$ and the specific forms of the functions are determined by the structures and rate constants of the M chemical reactions. As an example consider the set of reactions

$$PA$$
(6)
$$AB$$
(7)
$$A+2B3B$$
(8)
$$BC$$
(9)
The reactions represent the conversion of the precursor species $`P`$ into a final product $`C`$ via a sequence of four reactions involving two intermediates $`A`$ and $`B`$ . The third reaction is autocatalytic as $`B`$ catalyses its own production. The second reaction represents the uncatalysed conversion of $`A`$ to $`B`$ and the last reaction shows that the catalyst $`B`$ decays into the product $`C`$ . The equations are assumed to be irreversible. Also, the concentration of the reactant $`P`$ is assumed to be constant over a reasonable period of time. This is possible if the initial concentration of $`P`$ is large. Let $`p_0`$ (constant), $`a`$ and $`b`$ denote the concentrations of the molecular species $`P`$ , $`A`$ and $`B`$ . The decay product $`C`$ does not participate in any further reaction and so does not influence the chemical kinetics. The rates of the four successive chemical reactions are $`k_0p_0`$ , $`k_ua`$ , $`k_cab^2`$ and $`k_db`$ respectively where $`k_0`$ , $`k_u`$ , $`k_c`$ , $`k_d`$ are the rate constants. The equations governing the chemical kinetics are

$$\frac{da}{dt}=k_0p_0k_cab^2k_ua$$
(10)
$$\frac{db}{dt}=k_ua+k_cab^2k_db$$
(11)
In the general scheme of p.d.e.’s shown in Eq. (5), $`N`$ = 2, i.e., there are two molecular species $`A`$ and $`B`$ participating in $`M`$ = 4 chemical reactions. $`X_1=a`$ and $`X_2=b`$ are the concentrations of the moleculeas $`A`$ and $`B`$. The r.h.s.’s of Eqs. (10) and (11) are the nonlinear functions $`f_1(X_1,X_2)`$ and $`f_2(X_1,X_2),`$ the nonlinearity arising from the autocatalytic term $`k_cab^2`$ . The rate constants together with $`p_0`$ constitute the control parameters $`R`$.

In the general case, imagine an abstract $`N`$ \- dimensional state space with axes $`X_1,\mathrm{}.,X_N`$ . The state of the system at any instant of time, say $`t_0`$, is given by the $`N`$ -component state vector $`X(t_0)`$. In the state space, this state is represented by a single point. The time evolution of the system gives rise to a trajectory in the state space. The trajectory may end up at a fixed point $`X^{}`$. At this point, the rates of change of all the variables in the system are exactly zero, i.e., the l.h.s.’s of the $`N`$ equations in Eq.(5) are zero. The system is said to be in the steady state at the fixed point. At this point, the state of the system remains unchanged as a function of time. The only way of changing the state of the system is to apply perturbations to it. A fixed point is stable if small perturbations around the point eventually damp out. The stable fixed point acts as an attractor to the states in its vicinity. The corresponding region in the state space is called the basin of attraction. The nonlinear dynamics may give rise to more than one fixed point. If there are two stable fixed points, the system is bistable, i.e., two stable steady states are possible. One can similarly define multistability.

The other long-term possibilities for the trajectory in the state space are a limit cycle and a strange attractor. In the first case, the trajectory goes towards a closed loop and eventually circulates around it forever. In physical terms, this corresponds to stable oscillations in the system. The strange attractor is a set of states to which the trajectory is confined, never stopping or repeating. Such aperiodic motion is often indicative of chaos in the system. We now discuss the role of the control parameters R (Eqs. (4) and (5)) in the nonlinear dynamics of a system. By varying these parameters, one can bring about changes in the qualitative structure of the dynamics. Such changes are known as bifurcations. For example, as a parameter is changed, a steady state can become unstable and replaced by stable oscillations. A system with one stable steady state changes over to multistability, i.e., the system can exist in multiple steady states. To give a simple example of bifurcation, consider the rate equation

$$\frac{dx}{dt}=\mu xx^2$$
(12)
There are two fixed points of this equation: $`x^{}`$ = 0 and $`x^{}=\mu `$. To determine the stability of the fixed points, one undertakes what is known as the linear stability analysis. One determines the time evolution of a small perturbation $`\delta x(t)(=x(t)x^{})`$ around the fixed point. By substituting $`x(t)=x^{}+\delta x(t)`$ in Eq.(12) and ignoring terms of the order of $`(\delta x(t))^2`$ , one obtains $`\delta x(t)e^{\mu t}`$ when $`x^{}`$ = 0. The fixed point is stable if $`\mu <0`$ since $`\delta x(t)`$ reduces to zero during time evolution. The fixed point is unstable if $`\mu `$ \> 0 and $`\mu _c`$ = 0 is the bifurcation point. If $`x^{}=\mu `$, then $`\delta x(t)e^{\mu t}`$ . Hence the fixed point is unstable if $`\mu `$ \< 0 and stable for $`\mu `$ \> 0. Different types of bifurcation are possible a detailed discussion of which is given in standard textbooks and reviews  on nonlinear dynamics.

If there is more than one stable fixed point, a switch-like behaviour is possible. In the case of bistability, the system remains in one stable state until a sufficiently large perturbation drives the system to the other stable state. The system continues to remain in the latter state even after the perturbation is removed. The $`\lambda `$-phage lysis-lysogeny network offers an example of bistability . The lytic and the lysogenic states are the two possible steady states. A transition from the lysogenic to the lytic state occurs on irradiating with ultra-violet light. In a gene regulatory network, a negative (positive) feedback implies that a gene product inhibits (promotes) its own level of activity. To give an example, a protein which represses the transcription of its own gene operates through negative feedback. It has been found that negative (positive) feedback increases (decreases) stability in gene regulatory systems . Real life gene regulatory networks are often complex. Some of the examples are the $`\lambda `$-phage lysis-lysogeny circuit, the regulatory network for the activation of the tumour-suppressor protein p53  and the bacteriophage T7 (another lytic phage which infects E.coli ) network . Computational modelling studies of these networks have been undertaken with a view to explain experimental results. The quantitative agreement between theory and experiment is most often not good. The reasons are two fold: the complex nature of the networks and the difficulty in carrying out actual experiments on them. Computational as well as mathematical modelling of simpler networks is more extensive. Such networks incorporate the essential features of their more complex counterparts. The models seek to explain experimental results at a qualitative level. There are also abstract mathematical models of gene expression/regulation which highlight the general principles and their outcomes. There are already some good reviews and books on the computational and mathematical modelling of gene regulatory networks . For the purpose of this review, we pick on just one example, that of a synthetic gene regulatory network which illustrates the importance of nonlinearity in the dynamics of the network.

Gardner et al  have constructed and tested a synthetic, bistable gene regulatory network based on the predictions of a simple mathematical model. The network is called a genetic toggle switch and consists of two repressors (proteins) and two promoters. The enzyme RNA polymerase (RNAP) binds to the promoter region of a DNA sequence to initiate the process of transcription. The initial binding of RNAP to a promoter can be prevented by the binding of a regulatory protein to an overlapping segment of DNA, called operator. The gene expression is off in this case. Fig. 2 shows a simple sketch of the toggle network. The two promoters are designated as $`P_L`$ and $`Ptrc2`$ . $`P_L`$ drives the expression of the $`lacI`$ gene and $`Ptrc2`$ that of the $`cI`$ gene. The $`lacI`$ and $`cI`$ genes express the proteins of the same names. The proteins mutually inhibit the production of each other, hence the name repressor. The $`lacI`$ proteins form tetramers and the tetramer binds to operator sites adjacent to the $`Ptrc2`$ promoter, blocking the transcription of the $`cI`$ gene in the process. The $`cI`$ proteins, when produced, form dimers. The repressor dimer cooperatively binds to the operator sites in the vicinity of the $`P_L`$ promoter. As a result, transcription of the $`lacI`$ gene is not possible.

The nonlinear dynamics of the toggle network are governed by the following two equations:

$$\frac{dU}{dt}=\frac{\alpha _1}{1+V^\beta }U$$
(13)
$$\frac{dV}{dt}=\frac{\alpha _2}{1+U^\gamma }V$$
(14)
where $`U`$ and $`V`$ are the concentrations of $`lacI`$ and $`cI`$ proteins respectively, $`\alpha _1`$ and $`\alpha _2`$ are the effective rates of synthesis of $`lacI`$ and $`cI`$ proteins, $`\beta `$ is the cooperativity of repression of the $`P_L`$ promoter and $`\gamma `$ the same in the case of the $`Ptrc2`$ promoter. Fig. 3 reveals the origin of bistability in the system. The nullclines $`\frac{dU}{dt}`$ = 0 and $`\frac{dV}{dt}`$ = 0 intersect at three points. These are the fixed points (steady states) of the dynamics. Two of the fixed points are stable and the third unstable. The bistability occurs provided $`\beta ,\gamma `$ \> 1 (cooperative repression of transcription) and the rates of synthesis of the two repressors are balanced. If the rates are not balanced, the nullclines intersect at a single point giving rise to a single stable steady state (monostability).

In the region of bistability, the two stable steady states correspond to (1) State 1 (high $`V`$ / low $`U`$ ) and (2) State 2 ( low $`V`$ / high $`U`$ ) respectively. There are two basins of attraction, one above the separatrix and the other below it. In the $`log(\alpha _1)`$ vs. $`log(\alpha _2)`$ parameter space, bifurcation lines separate the monostable and bistable regions . The size of the bistable region decreases on reducing the cooperativity of repression ($`\beta `$ and $`\gamma `$ ). The parameters $`\alpha _1,\alpha _2,\beta `$ and $`\gamma `$ act as the control parameter R changing which a transition (bifurcation) between monostability and bistability occurs.

In the region of bistability, the toggle is flipped between the stable states (States 1 and 2) using transient chemical or thermal induction. The chemical agent isopropyl-$`\beta `$ -D-thiogalactopyranoside (IPTG) can bind to $`lacI`$ tetramers. As a result, the latter cannot bind to the operator region in the neighbourhood of the promoter $`Ptrc2`$, i.e., $`lacI`$ can no longer repress the production of the $`cI`$ proteins. Suppose the bistable system is originally in the stable State 2 (high $`U`$ ( $`lacI`$ ), low $`V`$ ( $`cI`$ )). On the induction of IPTG, the concentration of the $`cI`$ proteins increases as a function of time. The $`cI`$ proteins in their turn repress the production of $`lacI`$ proteins the concentration of which begins to fall. The dynamics ultimately leads the system to the other fixed point (State 1). The system remains in this stable steady state (low $`U`$ / high $`V`$ ) even after the removal of the IPTG stimulus. How can the toggle flip back to State 2 ? This is achieved by using a temperature-sensitive $`cI`$ protein in the network. The degradation rate of this protein increases as temperature is raised. On raising the temperature to $`42^{}C`$ (actual experiment), the concentration of $`cI`$ proteins starts to fall. Since repression is less, the concentration of $`lacI`$ proteins starts to go up.

The system finally reaches the fixed point corresponding to the stable steady State 2. After the steady state is reached, the temperature of the system is reduced ($`32^{}C`$ in the experiment). The system continues to remain in the steady State 2. A full cycle of the switching process is now completed. The actual construction of the toggle switch has been accomplished in E.coli using the standard tools of molecular biology . There is a reasonable agreement between the theoretical predictions based on Eqs.(13) and (14) and the results obtained from experiments on the synthetic toggle network. The design of the network relies significantly on theoretical inputs like identification of the region of bistability, increasing the cooperativity in repression ($`\beta `$ and $`\gamma `$ ) to achieve bistability over a wider region in parameter space etc. As a practical device, the toggle switch may have applications in biotechnology, biocomputing and gene therapy. As a cellular memory unit, the toggle provides the basis for “genetic applets” which are self-contained, programmable synthetic gene networks used in the control of cell functions. In parallel with the toggle work, another synthetic network, the repressilator has been designed and tested . The repressilator dynamics is again nonlinear and give rise to oscillations in the concentrations of the cellular proteins. The design of the network is based on a simple mathematical model of transcriptional regulation. The repressilator provides insight about the design principles of other oscillatory systems such as circadian clocks found in many organisms including cyanobacteria. The genetic toggle switch and the repressilator demonstrate that theoretical models can provide the design criteria for the actual construction of synthetic, gene regulatory networks. These simple networks have applications as practical devices and also help us to understand the functional properties of the more complex, naturally-occurring networks.

Nonlinear dynamics can give rise to various types of instability one of which is the Turing instability. In 1952, Turing  in a seminal paper proposed a mechanism for pattern formation in biological systems as well as the development of structure during the growth of an organism. Examples of biological patterns are the spots on the skin of a leopard, the stripes of a zebra, the arrangement of veins on the leaves of a tree etc. Structure formation is initiated by the process of cell differentiation, an example is the emergence of limbs in an organism during the growth of the organism from the featureless embryonic stage. The Turing mechanism involves both reaction as well as diffusion processes. To illustrate the mechanism, consider two chemical agents, the activator and the inhibitor. The activator is autocatalytic, i.e., it promotes its own production as well as that of the inhibitor. The inhibitor, as the name implies, is antagonistic to the activator and represses its production. Both the chemicals can diffuse but the inhibitor has a much larger diffusion coefficient. Consider a homogeneous distribution of the activator and the inhibitor in the system. Increase the concentration of the activator by a small amount in a local region. This gives rise to further increases in the local concentrations of the activator and the inhibitor. The inhibitor quickly diffuses to the surrounding region and prevents the activator from reaching there. Thus, in the steady state, islands of high activator concentration exist in a sea of high inhibitor concentration. The islands constitute what is known as the Turing pattern. Diffusion in general smooths out concentration differences in a system but the Turing process involving both reaction and diffusion gives rise to a steady pattern of concentration gradients. There is now increasing evidence that chemical gradients play a crucial role in the formation of patterns and cell differentiation in biological systems. To give an example, the protein bicoid has been found to have a graded concentration distribution in the Drosophila melanogaster embryo. It is responsible for the organization of the anterior half of the fly and has been fully characterised . Many reaction-diffusion (RD) models have been proposed based on the Turing mechanism and some of these can reproduce the patterns observed in nature . The basic scale of a pattern is larger than the size of an individual cell and so the RD processes involve more than one cell. Cells possibly choose developmental pathways depending upon their location in the concentration gradient. Position-dependent activation of genetic swtches in the cells may constitute an important step in both pattern and structure formation. Direct evidence for this in terms of a detailed characterization of the genes involved and an identification of the actual biochemical reactions occurring in the cells, is, however, yet to be obtained. Turing patterns have so far been experimentally observed in certain chemical RD systems in the laboratory  and also in some biological systems .

## 4 Effect of stochasticity

As already mentioned in the Introduction, stochastic fluctuations in the dynamics of the gene regulatory network lead to a probabilistic selection of developmental pathways. The $`\lambda `$-phage lysis-lysogeny network  was discussed as an example. Fig. 4 shows some of the key components of the network. The complexity of the full network is captured in Figure 1 of Ref. . We confine our attention to the simpler network. It consists of two $`\lambda `$ -phage genes $`cI`$ and $`cro`$. The corresponding promoters are $`P_{RM}`$ and $`P_R`$ respectively. Transcription of the gene $`cI`$ ( $`cro`$ ) expresses the regulatory protein $`\lambda `$ repressor ( $`Cro`$ ). Both the proteins are capable of binding to the operator regions $`O_R1`$ , $`O_R2`$ and $`O_R3`$. They act antagonistically to control promoter activity. Transcription of the $`cI`$ gene, initiated from the promoter $`P_{RM}`$ , takes place whenever there is no protein of either type binding to $`O_R3`$ . The $`\lambda `$ repressor molecule has a dumbell shape and there is a tendency for two such molecules to bind and form a dimer. The operator region $`O_R1`$ has the highest affinity for the binding of $`\lambda `$ repressor dimer. The binding increases the affinity of $`O_R2`$ for a second repressor dimer, i.e., a cooperative binding of dimers to the operator regions $`O_R1`$ and $`O_R2`$ takes place. The $`\lambda `$ repressor has both negative and positive control. If the $`\lambda `$ repressor is present at $`O_R2`$ , transcription of the $`cro`$ gene is not possible. This is because the repressor covers part of the $`DNA`$ that a $`RNAP`$ molecule must have access to in order to recognize the promoter $`P_R`$ , bind to it and initiate the transcription of the $`cro`$ gene. The same repressor at $`O_R2`$ exhibits positive control in helping a RNAP molecule to bind to the promoter $`P_{RM}`$ and begin transcription of the $`cI`$ gene. The increase in the transcription rate is approximately tenfold . If $`O_R2`$ is not occupied by the repressor, the transcription rate of the $`cI`$ gene is low. The reason for the dramatic increase in the transcription rate is the following. The presence of a repressor dimer bound to $`O_R2`$ leads to an increased affinity of $`P_{RM}`$ for $`RNAP`$ because the polymerase is held at $`P_{RM}`$ not only by its contacts with the $`DNA`$ but also due to the protein-protein contact with the repressor. In summary, a repressor dimer bound to $`O_R2`$ , represses transcription from $`P_R`$ but promotes transcription at $`P_{RM}`$ .

The $`cro`$ gene is transcribed only when the operator region $`O_R3`$ is either empty or has $`Cro`$ dimer bound to it. The transcription of the $`cI`$ gene cannot take place if the $`O_R1`$ and $`O_R2`$ regions are occupied by either protein, $`\lambda `$ repressor and $`Cro`$. In the lysogenic state, all the phage genes are off except for one gene $`cI`$ which produces the protein $`\lambda `$ repressor. The protein in turn binds to the operators $`O_R1`$ and $`O_R2`$ in the form of dimers and activates the transcription of its own gene at $`P_{RM}`$. The bound $`\lambda `$ repressor dimers further prevent transcription initiation at $`P_R`$ . Irradiation of the lysogen with ultra-violet light inactivates $`\lambda `$ repressor making the synthesis of the second regulatory protein $`Cro`$ possible. $`Cro`$ promotes lytic growth and competes with the $`\lambda `$ repressor in occupying the same operator sites. Increased $`Cro`$ production leads to a greater probability of $`Cro`$ binding at $`O_R3`$ which prevents the initiation of transcription at $`P_{RM}`$ . The concentration of the $`\lambda `$ repressor starts to fall down as a result. The concentration of $`Cro`$ proteins increases and when it reaches a level such that the operator regions $`O_R1`$ and $`O_R2`$ begin to be occupied, the transcription at $`P_R`$ is also halted. The switchover from the lysogenic to the lytic state is further possible by $`recA`$ -mediated degradation of the $`\lambda `$ repressor ( $`recA`$ is a catalytic protein ).

Arkin et al  have analysed the stochastic kinetics of the full $`\lambda `$-phage network which consists of more genes and regulatory elements than shown in Figure 4. Their detailed investigations show that fluctuations in the rates of gene expression give rise to random patterns of protein production in individual cells and wide diversity in instantaneous protein concentrations across cell populations. Each cell has two developmental pathways: lytic and lysogenic. The pathway selection depends upon which protein, $`\lambda `$ repressor or $`Cro`$, takes control of the operator region. If it is the $`\lambda `$ repressor, the lysogenic pathway is chosen. If the $`Cro`$ takes control, the lytic pathway is selected. Due to stochastic fluctuations, the concentrations of $`\lambda `$ repressor and $`Cro`$ vary considerably from cell to cell tipping the balance in favour of one or the other pathway. As a result, initially homogeneous cell populations can partition randomly into distinct lytic and lysogenic subpopulations. Arkin et al have constructed a stochastic kinetic model of the $`\lambda `$-phage circuit and based on model calculations predicted the fraction of infected cells selecting the lysogenic pathway at different phage:cell ratios. The theoretical results are consistent with the experimental results of Kourilsky . The kinetic model uses the stochastic formulation of chemical kinetics , stochastic mechanisms of gene expression  and a statistical-thermodynamical model of promoter regulation . Probabilistic selection of developmental pathways occurs in several other gene regulatory networks producing stochastic phenotypic outcomes. Some examples are given in Table 4 of Ref. .

We now describe the well-known Gillespie algorithm  which is incresingly being used by biologists in the stochastic kinetic approach to the study of gene expression and regulation in different systems Let us consider a system of $`N`$ chemicals participating in $`M`$ reactions $`R_\mu `$. The state of the system at any instant of time $`t`$ is represented as $`(X_1,\mathrm{},X_N)`$ where $`X_i`$ is the number of molecules of the ith chemical species. Two questions have to be answered to determine how the system evolves in time: (1) when will the next reaction occur and (2) what type of reaction will it be ? Let

$`C_\mu dt`$ = the probability that an $`R_\mu `$ $`(\mu =1,\mathrm{},M)`$ reaction occurs in the next infinitesimal time interval $`dt`$ for a particular combination of the reactant molecules. Let $`h_\mu `$ be the number of distinct combinations of molecules available in the state $`(X_1,\mathrm{}X_N)`$ for the $`R_\mu `$ reaction.
As an example, consider the reaction

$$A+BC$$
(15)
Let $`X_1`$ and $`X_2`$ be the number of molecules of types A and B respectively. Then $`h=X_1X_2`$ . Let

$`a_\mu dt=`$ $`h_\mu C_\mu dt`$ be the probability that an $`R_\mu `$ reaction occurs in time $`(t,t+dt)`$ given the system is in the state $`(X_1,\mathrm{},X_N)`$ at time $`t`$.
The reaction probability density function $`P(\tau ,\mu )d\tau `$ is the probability that given the state $`(X_1,\mathrm{},X_N)`$ at time $`t`$, the next reaction will occur in the infinitesimal time interval $`(t+\tau ,t+\tau +d\tau )`$ and will be an $`R_\mu `$ reaction,

$$P(\tau ,\mu )d\tau =P_0(\tau )a_\mu d\tau $$
(16)
where $`P_0(\tau )`$ is the probability that no reaction occurs in the time interval $`(t,t+\tau )`$ and $`a_\mu d\tau `$ is the subsequent probability that an $`R_\mu `$ reaction occurs in the time interval $`(t+\tau ,t+\tau +d\tau )`$. Now

$$P_0(\tau +d\tau )=P_0(\tau )\left[1\underset{\nu =1}{\overset{M}{}}a_\nu d\tau \right]$$
(17)
where the expression inside the bracket is the probability that no reaction occurs in time $`d\tau `$ from the state $`(X_1,\mathrm{}.,X_N)`$. Eq. (17) can be solved to obtain

$$P_0(\tau )=exp\left[\underset{\nu =1}{\overset{M}{}}a_\nu \tau \right]$$
(18)
Substituting for $`P_0(\tau )`$ in Eq. (16), one gets

$$P(\tau ,\mu )=a_\mu exp\left(a_0\tau \right)$$
(19)
if $`0\tau <\mathrm{}`$, $`\mu =1,\mathrm{},M`$ and $`P(\tau ,\mu )=0`$ otherwise.

$$a_\mu =h_\mu C_\mu ,(\mu =1,\mathrm{},M)$$
(20)
and

$$a_0=\underset{\nu =1}{\overset{M}{}}a_\nu $$
(21)
Now the goal is to generate a pair of random numbers $`(\tau ,\mu )`$ acording to the probability distribution (19). To do this, use the standard random number generator to obtain two random numbers from the uniform distribution in the unit interval. Take

$$\tau =\frac{1}{a_0}ln\left(\frac{1}{r_1}\right)$$
(22)
and $`\mu `$ is chosen to be the integer for which

$$\underset{\nu =1}{\overset{\mu 1}{}}a_\nu <r_2a_0\underset{\nu =1}{\overset{\mu }{}}a_\nu $$
(23)
The pair of numbers $`(\tau ,\mu )`$, (Eqs. (22) and (23)), belongs to the set of random pairs described by the probability density function $`P(\tau ,\mu )`$. For a rigorous proof of this see Refs. . Once $`(\tau ,\mu )`$ are known, put

$$t=t+\tau $$
(24)
and adjust the $`X_i`$ values according to the $`R_\mu `$ reaction. If the $`R_\mu `$ reaction is the one shown in Eq.(15), both $`X_1`$ and $`X_2`$ have to be decreased by 1 and $`X_3,`$ the number of molecules of C, increased by 1.

The input values at time $`t`$ = 0 are $`h_\nu ,C_\nu (\nu =1,\mathrm{}M)`$ and the initial values of $`X_i(i=1,\mathrm{},N)`$. The steps of the Gillespie algorithm are:

Step 1
Calculate $`a_\nu =h_\nu C_\nu (\nu =1,\mathrm{},M)`$ and $`a_0=_{\nu =1}^Ma_\nu `$.

Step 2
Generate $`r_1`$ and $`r_2`$ with the help of a uniform random number generator. Calculate $`\tau `$ and $`\mu `$ according to the formulae in Eqs. (22) and (23).

Step 3
Advance $`t`$ by $`\tau `$ (Eq.(24)) and adjust the $`X_i`$ values according to $`R_\mu `$. Then repeat the steps from Step 1 to further advance the system in time.

Recently, Kierzek et  al have used the Gillespie algorithm to study a stochastic kinetic model of prokaryotic gene expression. They explicitly considered ten biochemical reactions:

$$P+RNAPP_{}RNAP$$
(25)
$$P_{}RNAPP+RNAP$$
(26)
$$P_{}RNAPTrRNAP$$
(27)
$$TrRNAPRBS+P+EIRNAP$$
(28)
where $`P`$ denotes the promoter region of the gene and $`P_{}RNAP`$ the bound promoter-$`RNAP`$ complex. Reaction 2 (Eq.(26)) describes $`RNAP`$ dissociation and Reaction 3 the isomerization of “closed complex” to “open complex”, $`TrRNAP`$ is the activated $`RNAP`$promoter complex. Reaction 4 describes clearance of promoter region by $`RNAP,`$ $`EIRNAP`$ stands for $`RNAP`$ transcribing the gene and synthesizing the $`mRNA`$ molecule and $`RBS`$ is the ribosome binding site on $`mRNA`$. The other reactions are:

$$Ribosome+RBSRibRBS$$
(29)
$$RibRBSRBS+Ribosome$$
(30)
$$RibRBSEIRib+RBS$$
(31)
$$RBSdecay$$
(32)
$$EIRibprotein$$
(33)
$$Proteindecay$$
(34)
Reactions 5-10 (Eqs. (29)-(34)) describe translation, $`mRNA`$ decay and protein degradation. Reaction 5 describes $`Ribosome`$ binding to $`RBS`$ , the bound complex is designated as $`RibRBS`$ . Reaction 6 is the dissociation of the bound complex. Reaction 7 describes $`Ribosome`$ binding site clearance, $`EIRib`$ is the $`Ribosome`$ which translates the $`mRNA`$. Reaction 8 has degradation of $`RBS`$ by the enzyme $`RNAaseE`$. $`RNAaseE`$ and $`Ribosomes`$ are in competition to occupy $`RBS`$. If $`RNAaseE`$ binds first then it initiates the degradation of $`mRNA`$ but does not interfere with the movement of the already bound $`Ribosomes`$ engaged in the process of translation. Every $`Ribosome`$ which successfully binds to the $`RBS`$ completes translation of the protein. Reaction 9 corresponds to the completion of protein synthesis. Reaction 10 represents protein decay. The stochastic rate constants $`C_\mu ^{}s`$ of the different reactions, needed as inputs to the Gillespie algorithm, can be calculated from the more familiar chemical rate constants listed in Kierzek et al’s paper . For first order chemical reactions, the stochastic rate constant is equal to the rate constant of a chemical reaction. For second order reactions, the stochastic rate constant is equal to the rate constant divided by the volume of the system (in this case a cell).

In the last part of this Section, we describe a cooperative stochastic model of gene expression proposed by us . As already explained in the Introduction, the model has been constructed to explain the bimodal distribution in gene expression observed in recent experiments. The model describes the transcription of a single gene with one promoter region. There is one operaor region to which a regulatory protein $`R`$ can bind. This prevents the binding of a $`RNAP`$ to the promoter so that transcription of the gene cannot be initiated. There is a finite probability that the bound $`R`$ molecule dissociates from the operator at any instant of time. $`RNAP`$ molecule then has a certain probability of binding to the promoter and initiating transcription.

Each of the possibilities described above actually involves a series of physico-chemical processes, a detailed characterization of which is not required for the model of gene expression proposed by us. We represent a gene by a one-dimensional lattice of $`n+2`$ sites. The first two sites represent the operator and promoter respectively. The lattice is a coarse-grained description of an actual gene. In reality, the operator and promoter regions may extend over a certain number of base pairs in the DNA and they can be overlapping or not. In our model, they are represented as single sites. Each of the other sites in the lattice represents a finite number of base pairs in the DNA molecule.

The different physico-chemical processes are lumped together into a few simple events which are random in nature. This lumping together avoids unnecessary complexity that has no bearing on the basic nature of the process. The operator $`(O)`$ and the promoter $`(P)`$ together can be in four possible configurations: $`10,01,00`$ and $`11`$. The numbers “$`1`$” and “$`0`$” stand for “occupied” and “unoccupied”. The configuration $`ij`$ describes the occupation status of $`O`$ $`(i)`$ and $`P`$ $`(j)`$. For example, the configuration $`10`$ corresponds to $`O`$ being occupied by a $`R`$ molecule and $`P`$ being unoccupied. Similarly, in the configuration $`01`$, $`O`$ is unoccupied and $`P`$ is occupied by a $`RNAP`$ molecule. Binding of $`R`$ and $`RNAP`$ molecules are mutually exclusive so that the configuration $`11`$ is strictly prohibited. Given a $`00`$ configuration at time $`t`$, the transition probabilities to configurations $`10`$ and $`01`$ at time $`t+1`$ are $`p_1`$ and $`p_2`$ respectively. The probability of remaining in the configuration $`00`$ is $`1p_1p_2`$. A $`10`$ configuration at time $`t`$ goes to a $`00`$ configuration at time $`t+1`$ with probability $`p_3`$ and remains unchanged with probability $`1p_3`$. We have assumed all the probabilities to be time-independent. The $`RNAP`$ molecule once bound to the promoter initiates transcription in the next time step, i.e., the $`01`$ configuration makes a transition to a $`00`$ configuration with probability $`1`$. The motion of $`RNAP`$ is in the forward direction and the molecule covers a unit distance (the distance between two successive lattice sites) in each time step. Once the molecule reaches the last site of the lattice, the transcription ends and a $`mRNA`$ is synthesized.

The second major feature of our model is the cooperative binding of $`RNAP`$ to the promoter, when an adjacent $`RNAP`$ molecule is present. This implies that there is a higher probability of binding of $`RNAP`$ to the promoter in one time step if another $`RNAP`$ molecule is present at the site next to the promoter. In our model, the probability of cooperative binding is $`p_4`$ which is larger than $`p_2`$. The probabilities $`p_1`$ and $`1p_1p_2`$ are changed to new values $`p_5`$ and $`1p_4p_5`$ respectively. Degradation of $`mRNA`$ is taken into account by assuming the decay rate to be given by $`\mu N`$, where $`N`$ is the number of $`mRNAs`$ present at time $`t`$. The number of $`mRNAs`$ produced as a function of time is studied by Monte Carlo simulation. For the sake of simplicity, we have not tried to simulate protein levels or enzymatic products thereof, i.e., we study gene expression upto the level of transcription ($`mRNA`$ synthesis). Since the number of protein molecules and converted products should be proportional to the number of $`mRNA`$ molecules, no loss of generality is introduced by this simplification. The lattice consists of 52 sites $`(n=50)`$. Stochastic events are simulated with the help of a random number generator. The updating rule of our cellular automaton (CA) model is that in each time step $`t`$, the occupation status ( $`0`$ or $`1`$) of each site (except for the $`O`$ site) at time $`t1`$ is transferred to the nearest-neighbour site towards the right. If the last site is $`1`$ at $`t1`$, a $`mRNA`$ is synthesized at $`t`$ and the number of $`mRNAs`$ increases by one. In the same time step, the configuration $`ij`$ of $`OP`$ is determined with the probabilities already specified. Thus, in each time step, the $`RNAP`$ molecule, if present on the gene, moves forward by unit lattice distance (progression of transcription) followed by the updating of the $`OP`$ configuration. Figure 5 shows the concentration $`[mRNA]`$ of $`mRNA`$ molecules in the cell as a function of time for the parameter values $`p_1=0.5,p_2=0.5,p_3=0.3,p_4=0.85,p_5=0.05`$ and $`\mu =0.4`$. Note that an almost four-fold increase in the probability of $`RNAP`$ binding is assumed due to cooperativity. The stochastic nature of the gene expression is evident from the figure, with random intervals between the bursts of activity. One also notices the presence of several bursts of large size. It is important to emphasize that the frequency of transition between high and low expression levels is a function of the parameter values chosen and may be low for certain parameter values. For the probability values considered, the two predominantly favourable states are when the gene expression is off (state 1) and when a large amount of gene expression takes place (state 2). In the absence of $`RNAP`$ binding, state 1 has greater weight but with the chance binding of $`RNAP`$ to the promoter (probability $`p_2`$ for this is small), the weight shifts to state 2 until another stochastic event terminates cooperative binding and the gene reverts to state 1. The probability of obtaining a train of $`N`$ successive transcribing $`RNAP`$ molecules is $`p_2p_4^{N1}(1p_4)`$. This is the geometric distribution function and the mean and the variance of the distribution are given by $`\frac{p_2}{1p_4}`$ and $`\frac{p_2(1+p_4p_2)}{(1p_4)^2}`$ respectively. For the probability values already specified, the simulation has been repeated for an ensemble of 3000 cells. For each cell, the time evolution is upto 10000 time steps. Figure 6 shows the distribution of the number $`N(m)`$ of cells versus the fraction $`m`$ of the maximal number of $`mRNA`$ molecules produced after $`10000`$ time steps. Two distinct peaks are seen corresponding to zero and maximal gene expression. Such a bimodal distribution occurs over a range of parameter values.

Some theories have been proposed so far to explain the so-called “all or none” phenomenon in gene expression . These theories are mostly based on an autocatalytic feedback mechanism, synthesis of the gene product gives rise to the transport or production of an activator molecule. While such processes are certainly possible, the bimodal distribution is a much more general phenomenon and has now been found in many types of cells, from bacterial to eukaryotic and for different types of promoters . The two major features of the model of gene expression that we have proposed are stochasticity and cooperative binding of $`RNAP`$. There is by now enough experimental evidence of stochasticity in gene expression. Our suggestion of cooperative binding of $`RNAP`$ is novel and there is no direct experimental verification of the proposal as yet. There are some recent experiments which provide indirect evidence and these are discussed in Ref. .

## 5 Concluding remarks

In this review we have given an elementary introduction to some of the major aspects of biological networks, namely, topological characteristics, nonlinear dynamics and the role of stochasticity in gene expression and regulation. The main aim of the review is to highlight the usefulness of interdisciplinary approaches in the study of both natural and synthetic biological networks. Some important features of such networks have not been discussed in the review. One of these is the operational reliability of networks in spite of randomness in basic regulatory mechanisms. Many regulatory pathways do have highly predictable outcomes even when stochastic fluctuations are considerable. Cells adopt various strategies like populational transcriptional cooperation, checkpoints to ensure that cascaded events are appropriately synchronised, redundancy and feedback to achieve regulatory determinism. Some of these ideas are discussed in Ref. . The complexity of biological networks raises the question of their functional stability. In particular, the issue of interest is the sensitivity of networks to variations in their biochemical parameters. Barkai and Leibler  have studied a biochemical network responsible for bacterial chemotaxis and shown that the functional properties of the network are robust, i.e., relatively insensitive to changes in biochemical parameters like reaction rate constants and enzymatic concentrations. Bialek  has shown that extremely stable biochemical switches can be constructed from small numbers of molecules though intuitively one expects such systems to be prone to instability due to the inherent noise.

Metzler  in a recent paper has shown that spatial fluctuations in the distribution of regulatory molecules play a non-trivial role in genetic switching processes. Apart from internal stochastic fluctuations, external noise originating in random variations in the environment or in the externally set control parameters, may affect the functioning of a biological network. Hasty et al  have proposed a synthetic genetic network in which external noise is utilised to operate a protein switch (short noise pulses are used to turn protein production “on” and “off”). In another novel application, external noise is used to amplify gene expression, i.e., protein production by a considerable amount.

Genetic networks with many components are difficult to analyze using conventional techniques. Many parallels have been drawn in the functioning of genetic and electrical circuits . In electrical engineering, there are well developed techniques of circuit analysis which can be used to characterise the operation of complex electrical networks. Some of these techniques are increasingly being used to study genetic networks. Engineers are familiar with some of the design principles of biological networks. Rapid transitions between the two stable states of a system can be brought about by positive feedback loops. Negative feedback loops control the value of an output parameter to be within a narrow range even if there are wide fluctuations in the input. Coincidence detection systems activate an output provided two or more events occur simultaneously. Parallel connections enable a device to remain functional in the event of failures in one of the lines. One can give analogous examples from biology. One set of positive feedback loops is responsible for the rapid transition of cells into mitosis (division of cell nucleus), another set brings about the exit from mitosis in an irreversible manner. Gene transcription in eukaryotes involve coincidence detection. A $`mRNA`$ can be produced only if the promoters regulating gene expression are occupied by the different transcription factors. These examples indicate that general principles govern the functioning of genetic and electrical networks though there are other aspects of such networks which are not common to both. Biological networks constitute a field of research the interdisciplinary nature of which will become more evident as we progress into the twenty first century.
Acknowledgement: The Author thanks Subhasis Banerjee for help in drawing the figures.


format: text

text:
# 1 Introduction

## 1 Introduction

We will study the determination of a steeply falling spectrum when the energy measurement errors show long, asymmetrical, non-gaussian tails. We will see that for a steep spectrum the effect of long tail error functions can be dramatic: spectrum slopes, cutoff or other features could be changed or even erased, depending on the degree of deviation of the error function (energy resolution) from a gaussian distribution.

This problem, well known to nuclear spectroscopists working near the end point of a spectrum  is, however, often ignored in areas such as Cosmic Ray physics and we therefore choose to illustrate it by examining the spectrum of Ultra High Energy Cosmic Rays (UHECR) as determined by ground arrays, addressing the crucial question of the existence or absence of a cutoff in the spectrum, known as the Greisen-Zatsepin-Kuzmin (GZK)cutoff .

## 2 Spectrum reconstruction with non-gaussian errors

We will consider four cases for the reconstruction error function. These cases illustrate the different possibilities.

Let’s start with the simplest, well-known, case of a gaussian error reconstruction function. Assume that events of energy $`E^{}`$ are reconstructed with energy $`E`$ with a probability given by

$$P(E,E^{})=N\mathrm{exp}(\frac{(EE^{})^2}{2\sigma ^2}),$$
(1)
where $`\sigma `$ is constant. Then, the effect on an initial power law spectrum, $`\varphi =AE^\gamma `$, will be simply given by the convolution of $`\varphi `$ with the error function

$$\stackrel{~}{\varphi }(E)=_0^{\mathrm{}}𝑑E^{}P(E,E^{})\varphi (E^{}).$$
(2)
In this case the integral can be carried out analytically but is not very illuminating. For high energy $`E`$ it can be expanded to give

$$\stackrel{~}{\varphi }(E)=AE^\gamma \left(1+\gamma (\gamma 1)\frac{\sigma ^2}{2E^2}+\mathrm{}\right).$$
(3)
At large energies, the effect of the reconstruction energy on the reconstructed flux is very small. This is what is expected for gaussian errors and therefore one can safely neglect this effect. The same occurs for any other initial spectrum. The reconstructed spectrum will differ from the original one by a function which goes to zero rapidly with increasing energy.

Now let’s consider the log-gaussian distribution

$$P(E,E^{})=N\mathrm{exp}(\frac{\mathrm{log}(E/E^{})^2}{2\mathrm{\Delta }\xi ^2}),$$
(4)
where $`N`$ is a normalization constant and $`\mathrm{\Delta }\xi `$ is the standard deviation. Then the convolution of $`P`$ with the power like initial flux is

$$\stackrel{~}{\varphi }(E)=_0^{\mathrm{}}𝑑E^{}\varphi (E^{})P(E,E^{})=AE^\gamma \mathrm{exp}(\frac{\mathrm{\Delta }\xi ^2}{2}(\gamma ^21)),$$
(5)
where we are assuming that $`\mathrm{\Delta }\xi `$ is independent of energy. We see therefore that for a log-gaussian error with a constant standard deviation the effect of the convolution is to change the normalization. This is rather different from the result of Eq.3, where the change was negligible at large enough energies. Here the flux increases (for $`\gamma >1`$) by a constant factor. Notice also that the enhancement factor depends strongly on the initial spectrum index, $`\gamma `$. If we put $`\gamma =3`$ and $`\mathrm{\Delta }\xi 15`$% then $`\stackrel{~}{\varphi }/\varphi 1.09`$, for $`\mathrm{\Delta }\xi 30`$% then $`\stackrel{~}{\varphi }/\varphi 1.43`$. Alternatively, for $`\gamma =3`$, we get a factor two enhancement for $`\mathrm{\Delta }\xi 0.4`$.

In order to see what is the effect of very long tails on the reconstruction of events, let us consider an error function with a long tail

$$P(E,E^{})=N(\frac{E}{E^{}})^\alpha (1+\frac{E}{E^{}})^{\alpha \beta },$$
(6)
where $`N`$ is a normalization and $`\alpha `$ and $`\beta `$ are constants. The energy is reconstructed with a power law distribution of slope $`\alpha `$ and $`\beta `$ for small and large reconstructed energies. The convolution with the original spectrum gives as before

$$\stackrel{~}{\varphi }(E)=AK(\alpha ,\beta ,\gamma )E^\gamma ,$$
(7)
where $`K`$ is a constant which depends on the indexes $`\alpha ,\beta ,\gamma `$. As in the log-gaussian case the normalization changes but now the effect can be dramatic, depending on the values of $`\alpha `$ and $`\beta `$. For $`\alpha =4`$ and $`\beta =7`$ we find $`\stackrel{~}{\varphi }/\varphi 1.5`$. The log-gaussian distribution can be thought as the limit when $`\alpha `$ and $`\beta `$ are large.

In both cases the normalization is strongly affected. This is easily understood, the long tails on the error functions favor the more numerous low energy events to be reconstructed at higher energy, giving a higher flux. In all these cases we are assuming that there are no systematical energy shifts: the average reconstructed energy for a fixed initial energy is assumed to be equal to the initial energy. The existence of systematic errors in the reconstruction of energy would worsen considerably our results since in this case the change on the normalization of the spectrum will be linear in $`\mathrm{\Delta }\xi _{\mathrm{syst}}`$ rather than quadratic. Also we are assuming that the energy reconstruction function has no strong energy dependence.

Our last example of an error function with long tails is the Moyal distribution which constitutes an approximation to the Landau distribution describing the fluctuations in energy losses of an ionizing particle passing through matter .

$$P(E,E^{})=N\mathrm{exp}(\frac{1}{2}(\frac{(EE^{})}{\sigma }+e^{\frac{(EE^{})}{\sigma }})).$$
(8)
We will return to this distribution at the end of the paper when we look in detail at the problem of the GZK cutoff in the UHECR spectrum. Now let’s consider an initial spectrum with some feature. In order to keep the analysis simple, let’s assume a spectrum with a exponential cutoff.

$$\varphi (E)=AE^\gamma \mathrm{exp}(\frac{E}{E_c}),$$
(9)
where $`E_c`$ is the cutoff energy. The effect of a gaussian error would be again the original flux times a function going rapidly to 1 with energy, i.e. the original flux is not modified. However for a log-gaussian error or for a power law error we will have strong modification of the flux. In both cases at low energy ($`EE_c`$ ) we will have the previous result, the reconstructed flux is modified by a constant factor. At high energies, $`EE_c`$, the effect is different. In the log-gaussian error function we have

$$\stackrel{~}{\varphi }(E)=AE^\gamma \mathrm{exp}(\mathrm{\Delta }\xi ^2/2(\gamma ^21))e^{E/\stackrel{~}{E}_c},$$
(10)
where now $`\stackrel{~}{E}_c=E_c\mathrm{exp}(\gamma \mathrm{\Delta }\xi ^2)`$ is the new, smeared, cutoff energy. Again, this effect can be large due to the non linear dependence on $`\gamma `$ and $`\mathrm{\Delta }\xi `$; for $`\mathrm{\Delta }\xi =0.3`$ and $`\gamma =3`$ $`\stackrel{~}{E}_c=1.3E_c`$, for $`\mathrm{\Delta }\xi =0.4`$ $`\stackrel{~}{E}_c=1.62E_c`$.

In the case of a power law error function the effect is even more dramatic. At high energy, $`EE_c`$, the convolution of the error function with the exponential cutoff flux is

$$\stackrel{~}{\varphi }(E)=AE^\beta E_c^{\beta \gamma }K^{}(\alpha ,\beta ,\gamma ),$$
(11)
i.e. the cutoff is completely washed out and replaced by a spectrum with a constant slope. The new slope is solely dependent on the error function and not on the original spectrum.

We want now to argue that the non-gaussian error functions given above are realistic. There are several reasons why one expects long tails in the reconstruction of events in cosmic ray experiments. In ground array experiments the energy is estimated by measuring the density of particles in the shower. Usually parameters such as $`\rho (600)`$, the density of particles at $`600`$ m from the core, are used. It is found that the energy of the shower scales as $`E\rho (600)^\kappa `$ and $`\kappa 1`$. But this quantity has non-gaussian uncertainties. First, to calculate $`\rho (600)`$ (or any similar parameter) the core position and the arrival direction of the shower are needed. Given an arrival direction the $`\rho (600)`$ is corrected to ”zero degrees” arrival direction by an equation like

$$\rho (600)_0=\rho (600)\mathrm{exp}(s_0(\mathrm{sec}(\theta )1)),$$
(12)
where $`\theta `$ is the zenith angle and $`s_0`$ is a constant. Even if the arrival direction is reconstructed with a gaussian probability function the resulting $`\rho (600)_0`$ is not. However, this effect is, we believe, small at all but the largest zenith angles. The effect on the core reconstruction is more important. The density of particles at a distance $`r`$ from the shower core can be parameterized by the NKG lateral distribution function

$$\rho (r)=Kr^\eta (1+r/R)^\eta ^{},$$
(13)
where $`K,R,\eta ,\eta ^{}`$ are constants. Then, the error in the density due to the error in the determination of the core position would be generally non-gaussian.

Finally, any indirect measurement of cosmic rays based on shower development is subjected to shower fluctuations. It is well known that the total number of particles at fixed depth have large non-gaussian fluctuations, the total number of particles r.m.s scales as the energy. The same occurs for $`\rho (600)`$, since it is related to the total number of particles. In figure 1 we can see the distribution of the total number of particles at ground level in a shower of fixed energy. We can appreciate that the distribution of the number of particles is far from a gaussian and has long tails. In figure 2 we show the distribution of the density at 600 meters for showers of fixed energy. More than 20000 protons showers of energy $`10^{18}`$ were simulated with the Aires Monte Carlo code . As before the distribution is non-gaussian with tails extending to more than 1 order of magnitude above the average.

This was known for a long time but rarely used in the analysis of cosmic ray physics . Apparently, it is a generic phenomena in probability theory and does not contradict the central limit theorem .

## 3 The GZK Cutoff

Our results are not directly applicable to any real experiment. Rather, experiments should ascertain what are their error reconstruction functions and take them into account in the calculation of the spectrum. Also, events with large fluctuations could be cut off by other methods, improving therefore the energy reconstruction and avoiding these unwanted effects. Particularly important would be to get rid of any power like tail in the reconstruction error functions, since they have the most devastating effects on the spectrum reconstruction. On the other hand this effect could, at least in part, explain the current experimental situation where the normalization of the cosmic ray flux for different experiments is different and the existence of the GZK cutoff is controversial. Values of $`\mathrm{\Delta }\xi 3040`$ % are usual in cosmic ray physics and as said before with such values there are important effects on the shape and normalization of the resulting spectra.

The current experimental situation is aggravated by the low statistics in the high energy part of the spectrum. The Agasa experiment has measured $`58`$  events above $`4\times 10^{19}`$ eV, of them 8 are above $`10^{20}`$ eV. Using a realistic spectrum obtained from an uniform cosmological distribution of sources  and using a log-gaussian error function with $`\mathrm{\Delta }\xi =0.4`$ we obtain that for a total of 58 events above $`4\times 10^{19}`$ eV, the probability of having 8 or more events above $`10^{20}`$ eV is 2%. The average number of expected events would be 3. In figure 3 the convolution of different spectrum with energy error functions is shown. In the case of a $`\mathrm{\Delta }\xi =0.3`$ and the cosmological distribution of sources of Teshima and Yoshida  one can see that at $`E=200`$ EeV only a 0.1 reduction on the flux is expected.

The number of AGASA events above the GZK cutoff seems to establish its absence. However, if we introduce non-gaussian errors, we stress once again, the observation of such events becomes compatible with the GZK cutoff. To illustrate quantitatively the possibility of not seeing the GZK cutoff because of non-gaussian errors in the energy determination we consider the Moyal distribution of equation 4 (with approximately the same half-width of a gaussian distribution of a given variance, as seen in fig. 4) to each energy sampled from the AGASA spectrum with a exponential cutoff, that is a power law spectrum with an index of 2.78 multiplied by an exponentially falling spectrum starting at $`4\times 10^{19}`$ eV.

These errors alter the shape of the original spectrum and the expected number of events above the GZK cutoff . We summarize the results of this analysis in fig. 5.

We can compare the probability of finding events above GZK for measurements affected by gaussian errors (10%, 20%, and 30%) with measurements affected by Moyal-distributed errors (with the same half-width of gaussians with standard deviation of 10%, 20%, and 30% of energy). We see that considering only gaussian errors we cannot explain the observed number of events above GZK by AGASA (indicated by an arrow in figure 5) and would be led to conclude that there is a cutoff in the UHECR spectrum. However, for Moyal-distributed errors (30% of energy) we have an appreciable probability of finding more than 47 events above the GZK cutoff. We have shown in this paper that the experimental resolution curve when deviating from a gaussian distribution severely impacts the shape and normalization of a steeply falling spectra near its end point. Spectral features such as a cutoff may disappear from the reconstructed spectrum due to the effect of long tails in the resolution(error)function of the experiment.

Acknowledgments: We thank discussions with Máximo Ave and Enrique Zas. One of us (RAV) thanks the Instituto de Fisica Gleb Wataghim for kind invitation when this work took place and to FAPESP, Xunta de Galicia (XUGA-20604A98), and to the ”Ramón y Cajal” fellowship program for financial support. LGS and COE thank FAPESP for its financial support as well as the Brazilian Research Council-CNPq (COE).


text:
# Elementary vortex pinning potential in a chiral 𝑝-wave superconductor

## Abstract

The elementary vortex pinning potential is studied in a chiral $`p`$-wave superconductor with a pairing $`𝐝=\overline{𝐳}(\overline{k}_x\pm i\overline{k}_y)`$ on the basis of the quasiclassical theory of superconductivity. An analytical investigation and numerical results are presented to show that the vortex pinning potential is dependent on whether the vorticity and chirality are parallel or antiparallel. Mutual cancellation of the vorticity and chirality around a vortex is physically crucial to the effect of the pinning center inside the vortex core.

Much attention has been focused on the vortex pinning in type-II superconductors. The vortex pinning plays an important role on various vortex-related quantities and phenomena such as the critical current and the hysteresis of the magnetization in superconductors under magnetic fields. The characteristics of the vortex-related phenomena are of particular interest in unconventional superconductors with multiple components of the superconducting order parameter. In such superconductors, multiple states of superconducting order can coexist. Accordingly there appear multiple kinds of vortex structure, where the nature of the vortex pinning can be dependent on the microscopics of the superconducting order.

One of the superconductivity with multiple components of the order parameter is the chiral $`p`$-wave one, $`𝐝=\overline{𝐳}(\overline{k}_x\pm i\overline{k}_y)`$, which is composed of two degenerate pairing states $`\overline{k}_x`$ and $`\overline{k}_y`$ and breaks the time-reversal symmetry. This superfluid 3He-A type of chiral $`p`$-wave pairing state has been anticipated in a layered ruthenate superconductor Sr2RuO4. While the identification of the genuine superconducting pairing of this material is still open to further discussion, that chiral $`p`$-wave pairing has the simplest and essential form and has attracted a great deal of attention. The vortices for Sr2RuO4 have been investigated intensively. In the context of the vortex pinning, we will see a rich physics contained in that chiral $`p`$-wave pairing state.

In this paper, we investigate the elementary vortex pinning potential in the chiral $`p`$-wave superconductor with the pairing $`𝐝=\overline{𝐳}(\overline{k}_x\pm i\overline{k}_y)`$. A point-like pinning center and a single vortex with vorticity perpendicular to conduction layers in a layered superconductor are considered. We show that the vortex pinning potential depends on the sense of the chirality of the Cooper pairs relative to the vorticity of the vortex in the chiral $`p`$-wave superconductor. First we analytically discuss the interplay between the chirality and vorticity to explain the mechanism of the chirality dependence of the vortex pinning potential. We then present numerical result for the vortex pinning potential obtained from self-consistent order parameters. Our numerical result confirms the analytical one. The chirality dependence of the vortex pinning would have influence on the hysteresis of the magnetization and the distribution of the magnetic field in samples, which might be observed with SQUID and magneto-optical imaging techniques.

To investigate the vortex pinning, we use the quasiclassical theory of superconductivity. We start with the Eilenberger equation for the quasiclassical Green function in the absence of the pinning,

$$\widehat{g}_{\mathrm{imt}}(i\omega _n,𝐫,\overline{𝐤})=i\pi \left(\begin{array}{cc}g_{\mathrm{imt}}& if_{\mathrm{imt}}\\ if_{\mathrm{imt}}^{}& g_{\mathrm{imt}}\end{array}\right),$$
(1)
namely,

$`iv_\mathrm{F}\overline{𝐤}\widehat{g}_{\mathrm{imt}}+[i\omega _n\widehat{\tau }_z\widehat{\mathrm{\Delta }},\widehat{g}_{\mathrm{imt}}]=0,`$ (2)
where the order parameter is $`\widehat{\mathrm{\Delta }}(𝐫,\overline{𝐤})=\left[(\widehat{\tau }_x+i\widehat{\tau }_y)\mathrm{\Delta }(𝐫,\overline{𝐤})(\widehat{\tau }_xi\widehat{\tau }_y)\mathrm{\Delta }^{}(𝐫,\overline{𝐤})\right]/2`$ and $`\widehat{\tau }_i`$ the Pauli matrices. The Eilenberger equation (2) is supplemented by the normalization condition $`\widehat{g}_{\mathrm{imt}}(i\omega _n,𝐫,\overline{𝐤})^2=\pi ^2\widehat{1}`$, and the commutator is $`[\widehat{a},\widehat{b}]=\widehat{a}\widehat{b}\widehat{b}\widehat{a}`$. The vector $`𝐫=(r\mathrm{cos}\varphi ,r\mathrm{sin}\varphi )`$ is the center of mass coordinate, and the unit vector $`\overline{𝐤}=(\mathrm{cos}\theta ,\mathrm{sin}\theta )`$ represents the relative coordinate of the Cooper pair. The cylindrical Fermi surface is assumed. We use units in which $`\mathrm{}=k_\mathrm{B}=1`$.

Following Thuneberg et al., the effect of the pinning is introduced to the quasiclassical theory of superconductivity as follows. The quasiclassical Green function $`\widehat{g}`$ in the presence of a point-like non-magnetic defect situated at $`𝐫=𝐑`$ is obtained from the Eilenberger equation

$`iv_\mathrm{F}\overline{𝐤}\widehat{g}+[i\omega _n\widehat{\tau }_z\widehat{\mathrm{\Delta }},\widehat{g}]=[\widehat{t},\widehat{g}_{\mathrm{imt}}]\delta (𝐫^{}),`$ (3)
and the $`t`$ matrix due to the defect

$`\widehat{t}(i\omega _n,𝐫^{})={\displaystyle \frac{v}{D}}\left[\widehat{1}+N_0v\widehat{g}_{\mathrm{imt}}(i\omega _n,𝐫^{},\overline{𝐤})_\theta \right],`$ (4)
where $`𝐫^{}=𝐫𝐑`$, the denominator $`D=1+(\pi N_0v)^2\left[g_{\mathrm{imt}}_\theta ^2+f_{\mathrm{imt}}_\theta f_{\mathrm{imt}}^{}_\theta \right]`$, the average over the Fermi surface $`\mathrm{}_\theta =\mathrm{}𝑑\theta /2\pi `$, the normal-state density of states on the Fermi surface $`N_0`$, and we assume the $`s`$-wave scattering $`v`$ when obtaining Eq. (4). We define a parameter $`\sigma =(\pi N_0v)^2/\left[1+(\pi N_0v)^2\right]`$, which measures how strong the scattering potential of the defect is.

The free energy in the presence of the defect is, at the temperature $`T`$, given as

$`\delta \mathrm{\Omega }(𝐑)=N_0T{\displaystyle _0^1}𝑑\lambda {\displaystyle \underset{\omega _n}{}}{\displaystyle 𝑑\overline{𝐤}𝑑𝐫\mathrm{Tr}\left[\delta \widehat{g}_\lambda \widehat{\mathrm{\Delta }}_b\right]},`$ (5)
where $`\delta \widehat{g}_\lambda =\widehat{g}\widehat{g}_{\mathrm{imt}}`$ is evaluated at $`\widehat{\mathrm{\Delta }}=\lambda \widehat{\mathrm{\Delta }}_b`$, and $`\widehat{\mathrm{\Delta }}_b`$ is the order parameter in the absence of the defect. Equation (5) represents the difference in the free energy between the states with and without the defect, and then gives the vortex pinning potential $`\delta \mathrm{\Omega }(𝐑)`$.

For the chiral $`p`$-wave pairing state $`𝐝=\overline{𝐳}(\overline{k}_x+i\overline{k}_y)=\overline{𝐳}\mathrm{exp}(i\theta )`$, it is known that the order parameter around a single vortex, $`\mathrm{\Delta }_b(𝐫,\overline{𝐤})`$ $`\left[\mathrm{\Delta }_b(r,\varphi ;\theta )\right]`$, has two possible forms depending on whether the chirality and vorticity are parallel or antiparallel each other. One form is

$`\mathrm{\Delta }_b^+(r,\varphi ;\theta )=\mathrm{\Delta }_+(r)e^{i(\theta \varphi )}+\mathrm{\Delta }_{}(r)e^{i(\theta +\varphi )},`$ (6)
where the chirality and vorticity are antiparallel (Case I). The other is

$`\mathrm{\Delta }_b^{++}(r,\varphi ;\theta )=\mathrm{\Delta }_+(r)e^{i(\theta +\varphi )}+\mathrm{\Delta }_{}(r)e^{i(\theta +3\varphi )},`$ (7)
where the chirality and vorticity are parallel (Case II). Here, the vortex center is situated at $`𝐫=0`$, the dominant component $`\mathrm{\Delta }_+(r\mathrm{})=\mathrm{\Delta }_{\mathrm{BCS}}(T)`$, and the induced one $`\mathrm{\Delta }_{}(r\mathrm{})=0`$. Because of axisymmetry of the system, we can take $`\mathrm{\Delta }_\pm (r)`$ to be real.

First we analytically investigate the vortex pinning potential. We discuss the quantity $`\delta \mathrm{\Omega }(R=0)`$, where both the defect and the vortex center are situated just at the origin $`𝐫=0`$ ($`R|𝐑|`$). From the quasiclassical viewpoint, the quasiparticles inside the vortex core, subject to the Andreev reflection, run along straight lines called as quasiparticle paths. We consider the quasiparticle paths which go through the origin $`𝐫=0`$. On those paths with zero impact parameter, the position vector is parallel to the direction of the quasiparticle path (i.e., $`𝐫\overline{𝐤}`$), and therefore $`\varphi =\theta ,\theta +\pi `$. In this situation ($`\varphi =\theta `$), from Eqs. (6) and (7), the order parameter on the path is

$`\mathrm{\Delta }_b^+(r,\varphi =\theta ;\theta )=\mathrm{\Delta }_+(r)+\mathrm{\Delta }_{}(r)`$ (8)
in case I, and

$`\mathrm{\Delta }_b^{++}(r,\varphi =\theta ;\theta )=[\mathrm{\Delta }_+(r)+\mathrm{\Delta }_{}(r)]e^{2i\theta }`$ (9)
in case II. The cancellation between the chirality and vorticity occurs in Eq. (8) and not in Eq. (9). Of importance is the resultant difference in the phase factor of these order parameters.

On the basis of an analysis of the so-called zero-core vortex model in Ref. , the matrix elements of $`\widehat{g}_{\mathrm{imt}}`$ at the vortex center are approximately obtained as

$$g_{\mathrm{imt}}=\sqrt{\omega _n^2+|\stackrel{~}{\mathrm{\Delta }}|^2}\omega _n^1,f_{\mathrm{imt}}=\stackrel{~}{\mathrm{\Delta }}\omega _n^1,f_{\mathrm{imt}}^{}=\stackrel{~}{\mathrm{\Delta }}^{}\omega _n^1,$$
(10)
where $`\stackrel{~}{\mathrm{\Delta }}=\mathrm{\Delta }_b^{+\pm }(r\mathrm{},\varphi =\theta ;\theta )`$. Here, Eq. (10) is obtained with assuming that the amplitude of the order parameter is constant (i.e., zero core) around the vortex, which is the only approximation in this analysis. Inserting the order parameter of Eq. (8) into Eq. (10), we obtain the anomalous Green functions integrated over the Fermi surface as, in case I,

$$f_{\mathrm{imt}}_\theta =f_{\mathrm{imt}},f_{\mathrm{imt}}^{}_\theta =f_{\mathrm{imt}}^{}$$
(11)
because of the absence of any phase factors in Eq. (8), i.e., because of the cancellation between the chirality factor $`\mathrm{exp}(i\theta )`$ and the vorticity factor $`\mathrm{exp}(i\varphi )`$ in Eq. (6). On the other hand, in case II,

$$f_{\mathrm{imt}}_\theta =0,f_{\mathrm{imt}}^{}_\theta =0$$
(12)
because of the phase factor $`\mathrm{exp}(2i\theta )`$ contained in Eq. (9). The diagonal component of $`\widehat{g}_{\mathrm{imt}}_\theta `$ is $`g_{\mathrm{imt}}_\theta =g_{\mathrm{imt}}`$ both in cases I and II. Consequently, in case I, $`\widehat{g}_{\mathrm{imt}}_\theta =\widehat{g}_{\mathrm{imt}}`$ and we obtain $`[\widehat{t},\widehat{g}_{\mathrm{imt}}]=0`$ from Eq. (4). In case II, $`\widehat{g}_{\mathrm{imt}}_\theta \widehat{g}_{\mathrm{imt}}`$ and $`[\widehat{t},\widehat{g}_{\mathrm{imt}}]0`$ generally.

When $`[\widehat{t},\widehat{g}_{\mathrm{imt}}]=0`$, the Eilenberger equation (3) in the presence of the defect is identical to Eq. (2) (the equation in the absence of the defect), namely, the defect has no influence on the Green function and the free energy. From this and the above results of the analysis of the factor $`[\widehat{t},\widehat{g}_{\mathrm{imt}}]`$, we find that $`\delta \mathrm{\Omega }(0)=0`$ in case I when the chirality is antiparallel to the vorticity, and $`\delta \mathrm{\Omega }(0)0`$ in case II when the sense of the chirality is the same as that of the vorticity. It means that the vortex pinning depends on the chirality in the chiral $`p`$-wave superconductor.

The above analytical result is based on the zero-core vortex model, i.e., on the non-self-consistent (constant) amplitude of the order parameter. We next investigate the vortex pinning potential $`\delta \mathrm{\Omega }(R)`$ numerically with the self-consistent order parameters around the vortex which have the forms of Eqs. (6) and (7). As self-consistent amplitude $`\mathrm{\Delta }_\pm (r)`$ in Eqs. (6) and (7), we adopt numerical data which we have obtained in Ref.  by solving self-consistently the Eilenberger equation.

In Fig. 1, we show the numerical results for $`\delta \mathrm{\Omega }(R)`$ in the Born limit ($`\sigma 1`$) and the unitary limit ($`\sigma 1`$). We present those results for the chiral $`p`$-wave pairing and the isotropic $`s`$-wave one. As noted in Fig. 1, in the case of the $`s`$-wave pairing (dot-dashed lines), the difference in the free energy between the states with and without the defect, $`\delta \mathrm{\Omega }(R)`$, is equal to zero at $`R\mathrm{}`$ ($`R`$ is the distance between the vortex center and the defect). This is because the Anderson’s theorem  is valid far away from the vortex core. On the other hand, in the chiral $`p`$-wave pairing cases (solid and dashed lines), $`\delta \mathrm{\Omega }(R\mathrm{})`$ is finite and positive as seen in Fig. 1. The quantity $`\delta \mathrm{\Omega }(R\mathrm{})`$ is equal to the loss of the condensation energy due to the pair breaking effect of the defect far away from the vortex core (i.e., the breakdown of the Anderson’s theorem). As noted in Figs. 1(a) and 1(b), at $`T=0.8T_\mathrm{c}`$ (high temperature), the condensation energy loss in bulk $`\delta \mathrm{\Omega }(R\mathrm{})`$ dominantly contributes to the depth of the vortex pinning potential $`\delta \mathrm{\Omega }(R)`$, i.e., to the vortex pinning energy. From Figs. 1(a) and 1(b) it is noticed that the vortex pinning energies of the chiral $`p`$-wave pairing cases at a high temperature are about 10 times larger than those of the $`s`$-wave pairing case. This enhancement of the pinning effect is due to the breakdown of the Anderson’s theorem, and then it must be a common feature of unconventional superconductors. For example, in the case of high-$`T_\mathrm{c}`$ cuprates, this may be one of the reasons why small point defects such as Zn atoms and oxygen vacancies are efficient pinning centers.

As noted in Figs. 1(c) and 1(d), at $`T=0.2T_\mathrm{c}`$ (low temperature), the contribution of the vortex core ($`R0`$) to the depth of $`\delta \mathrm{\Omega }(R)`$ is nonzero in case II (dashed lines). Here, the contribution of the vortex core means the energy gain due to the presence of the scattering center in the vortex core. In contrast, the depth of $`\delta \mathrm{\Omega }(R)`$, i.e., the vortex pinning energy, is determined in case I (solid lines) only by the loss of the condensation energy far away from the vortex core. It is noticeable that certainly $`\delta \mathrm{\Omega }(R=0)`$ equals to zero in case I. This numerical result confirms the analytical one discussed above. The vortex pinning energy depends on whether the chirality and vorticity are antiparallel (solid lines) or parallel (dashed lines). Especially in the Born limit, the difference in the vortex pinning energy is eminent as noticed in Fig. 1(c), because in this limit the loss of the condensation energy in bulk is relatively small compared to the contribution of the vortex core to the depth of $`\delta \mathrm{\Omega }(R)`$.

In general, the two chiral states of cases I and II can coexist as domain structure in samples under magnetic fields. The spatial gradient of the magnetic field in a sample is proportional to the local strength of the vortex pinning in the critical state. In terms of the present chirality-dependent vortex pinning, the gradient inside the domain of the case-II state is predicted to be steeper than that inside the domain of the case-I state. This may be experimentally observed as a signature of the chiral state. Also the domain structure of the two chiral states depends on the hysteresis of applied magnetic field, and therefore the present chirality-dependent vortex pinning may affect the hysteresis curve during multiple cycles of the magnetization as observed in UPt3.

In the case of the usual winding-1 vortex $`\mathrm{\Delta }\mathrm{exp}(i\varphi )`$, the chiral “$`p`$-wave” pairing $`\overline{k}_x\pm i\overline{k}_y=\mathrm{exp}(\pm i\theta )`$ is essential for the cancellation between the chirality and vorticity. If winding-2 vortices $`\mathrm{\Delta }\mathrm{exp}(2i\varphi )`$ are realized in a chiral $`d`$-wave state $`\overline{k}_x^2\overline{k}_y^2\pm i\overline{k}_x\overline{k}_y=\mathrm{exp}(\pm 2i\theta )`$, the same kind of cancellation occurs.

We comment on the relation of the present vortex pinning to the superconducting gap structure in Sr2RuO4. In this material, it has been pointed out from experiments that the gap has line nodes and little in-plane anisotropy. Models for the gap structure consistent with those experimental facts were proposed, in which there existed horizontal line nodes perpendicular to the axis of the cylindrical Fermi surface. Now, for the present theory of the chirality-dependent vortex pinning, what is important is that the Fermi surface averages of the anomalous Green functions (i.e., the average of the order parameter except for the chiral part) are finite as in Eq. (11). The present chirality dependence of the vortex pinning does not occur, if the order parameters have sign changes on all Fermi surfaces relevant to superconductivity as $`\mathrm{\Delta }(𝐤)\mathrm{exp}(\pm i\theta )\mathrm{cos}(ck_z)`$. It occurs, if there are no sign changes as $`\mathrm{\Delta }(𝐤)\mathrm{exp}(\pm i\theta )|\mathrm{cos}(ck_z)|`$. In another case, the chirality dependence is expected to occur when the order parameter is nodeless on the major Fermi surface with dominant density of states, even if there are gap nodes and sign changes on the other minor Fermi surfaces.

In conclusion, we investigated the elementary vortex pinning potential $`\delta \mathrm{\Omega }(R)`$ on the basis of the quasiclassical theory of superconductivity. In the chiral $`p`$-wave pairing state, $`\delta \mathrm{\Omega }(R)`$ was dependent on the sense of the chirality relative to the vorticity at a low temperature. In terms of the present chirality-dependent vortex pinning, a theoretical analysis for anomalies in the hysteresis of the magnetization observed experimentally in Sr2RuO4 would be interesting and is left for future work.

We thank T. Tamegai, P. H. Kes, T. Kita, K. Kinoshita, and A. Maeda for helpful discussions, and J. R. Clem, E. V. Thuneberg, and M. L. Kulić for considerate correspondence. One of the authors (N.H.) also thanks M. Ichioka, K. Machida, M. Takigawa, N. Nakai, M. Matsumoto, S. H. Pan, and M. Fogelström for useful discussions and comments. Y.K. thanks M. Sigrist for useful discussions. This work is partly supported by Grant-in-Aid for Scientific Research on Priority Areas (A) of “Novel Quantum Phenomena in Transition Metal Oxides” (12046225) from the Ministry of Education, Science, Sports and Culture and Grant-in-Aid for Encouragement of Young Scientists from Japan Society for the Promotion of Science (12740203).