text:
# Anomalous relaxation kinetics of biological lattice-ligand binding models

## I Introduction

How does a system evolve towards its steady state? Sometimes the answer is quite simple and the relaxation process is merely an exponential decay. If the deviations from equilibrium are small Onsager’s regression hypothesis onsager31  asserts that the relaxation is governed by the same laws as the fluctuations in equilibrium. This hypothesis certainly fails for systems with an absorbing steady state such as simple models for diffusion-limited chemical reactions torney83 ; toussaint83 ; racz85 ; peliti86 ; nielaba92 . Here there are no fluctuations in the steady state but the approach towards the absorbing state is critical in the sense that it exhibits slow power-law decay and universal scaling behavior cardyreview . Most of these models are chosen to be mathematically transparent hoping that they still resemble some of the essential features of actual systems occurring in nature. Unfortunately, experimentally accessible systems where the above theoretical ideas can be tested explicitly have remained rare to date.

In this contribution we discuss the kinetics of some macromolecular assembly processes relevant for the formation of functional structures in cells. In particular, we are interested in the dynamics of ligand-substrate binding, where the substrate is a one- or two-dimensional lattice and the ligands are dimers or oligomers. Examples for such systems are illustrated in Fig. 1. Binding of dimeric myosin to actin filaments conibear98  can behave in a similar way as dimeric kinesin on microtubules. Similar interactions of large supramolecular biological polymers with protein ligands such as DNA with proteins or viruses with antibodies are also a quite intensive area of research. As will become clear in later chapters the kinetics of these systems shows interesting anomalous dynamics which is closely related to the mathematical models of simple diffusion-reaction systems discussed above.

Fig. 1a shows a schematic representation of a “decoration experiment” where dimeric motor enzymes (ligands) are deposited on their corresponding molecular track (substrate) mandelkow\_hoenger99 ; hoenger00 . For kinesin motors, these tracks are microtubules, hollow cylinders usually consisting of $`13`$ protofilaments, linear polymers composed of alternating $`\alpha `$\- and $`\beta `$-tubulin subunits. The kinesin binding sites are located on the $`\beta `$ subunits (dark spheres) which form a helical (wound-up rhombic) lattice with a longitudinal periodicity of $`8\mathrm{nm}`$. Kinesin is a mechanochemical enzyme, which transforms (through an isothermal stochastic process) the chemical energy obtained from the hydrolysis of adenosine triphosphate (ATP) into motion along microtubules. Roughly speaking, kinesin has the following building plan: two globular heads (also called motor domains) with the functionality of ”legs” are joined together by a coiled-coil region into a tail which can bind to some cargo. The typical size of these proteins is in the range of several tens of nanometers.

Decoration techniques are usually performed in the absence of ATP hydrolysis. Then the motor enzymes are “passive” ligands which bind and unbind from their tubulin binding sites but do not actively move along microtubules. These systems have traditionally been used in biophysical chemistry to investigate the structure and the binding properties of kinesin mandelkow\_hoenger99 ; hoenger00 , i.e. after waiting for the system to equilibrate the binding stoichiometry is measured and the structure is determined by cryo-electron microscopy followed by 3D image reconstruction. These investigations are key elements in understanding the binding patterns of kinesin motor domains under changing nucleotide conditions to get a complete picture on the different conformational states, which are involved in the kinesin walking mechanism schief2001 .

In the following sections we present a theoretical analysis of dimer adsorption-desorption kinetics with competing single and double bound dimers. Such analysis is extremely important for a quantitative analysis of decoration data vilfan2001a ; a theoretical analysis gives the binding stoichiometry in the equilibrium state in terms of binding constants for the first and second head of the dimer molecule. The dynamics of the approach to equilibrium is useful to estimate when an experimental system can be considered as equilibrated. Even more importantly, time-resolved decoration experiments (e.g. by using motor enzymes labelled by some fluorescent marker) combined with our theoretical analysis could provide new information about reaction rates which are to date not known completely. Understanding the kinetics of passive motors is undoubtedly a necessary prerequisite for studying the more complicated case of active motors at high densities capello\_preprint . The model is also interesting in its own right since it contains some novel features of non-equilibrium dynamics of dimer adsorption-desorption models.

Fig. 1b shows a schematic representation of tropomyosin binding to an actin filament vilfan2001c . Actin filaments are one of the major compontents of the cytoskeleton. Like microtubules, they contribute to the mechanical stability of the cell and serve as tracks for molecular motors from the myosin family. This function is especially pronounced in muscle cells where both actin and myosin form filaments that can slide between each other and thereby cause the muscle to contract. Tropomyosin plays the key role in the control of skeletal muscle contraction. It binds to actin along its binding sites for myosin motors. When calcium ions are released as a response to a nerve signal, they cause tropomyosin to shift laterally thereby clearing the binding sites and allowing myosin to bind and produce force. Tropomyosin binding to action has several features which make it an interesting model system of statistical physics; each molecule covers seven actin monomers and interacts strongly with other tropomyosin molecules. As a consequence gaps between bound molecules can take a long time to heal. Although the gap dynamics has not yet been measured experimentally, the amount of data gathered in other kinetic studies provides plenty of information on the model parameters and allows to make quantitative predictions about the relaxation time. As it turns out vilfan2001c , these relaxation times can be as large as hours, which is short enough compared to the lifetime of actin filaments in a muscle cell, but essential when planning experiments with actin and tropomyosin assembled in vitro, especially when studying the myosin regulation fraser95 .

The outline of this article is as follows. In Section II we define the model for decoration of microtubules with dimeric motors which can bind either with one or with two heads as first introduced in Ref. vilfan2001a . In Section III we determine analytically the equilibrium state of this model. We study the dynamics of the model in Section IV, which extends the results presented in Ref. vilfan2001b . Section V studies the dynamics of the two-dimensional model with interactions and Section VI the dynamics of the one-dimensional $`k`$-mer model, which is a more generalized version of the results described in Ref. vilfan2001c . In the final section we give a summary and an outlook on future challenges in the field.

## II Definition of the microtubule decoration model

Our model describes the experimental situation as found in most decoration assays. It starts with an empty tubulin sheet surrounded by a solution of double-headed kinesin molecules. The kinesin dimers can either attach with one or two heads onto binding sites located on $`\beta `$-tubulin. There seems to be convincing evidence that kinesin heads can bind on two adjacent binding sites only in longitudinal but not in lateral direction mandelkow\_hoenger99 ; hoenger00 . This introduces a strong uniaxial anisotropy and distinguishes the adsorption process of protein dimers from simple inorganic dimers. The attached heads can also detach at some rate. A schematic representation of the decoration process onto tubulin sheets is given in Fig. 2.

To begin with, let us take into account only steric (hard-core) interactions and for now neglect nearest neighbour attractive interactions. Then we are left with a one-dimensional problem of kinesin dimers decorating a single protofilament (one-dimensional lattice) as shown in Fig. 2 and defined as follows.111A Java applet with the simulation of our model can be found at http://www.ph.tum.de/~avilfan/relax/ . Kinesin is considered as a dimeric structure with its two heads tethered together by some flexible joint. Hence each dimer (kinesin protein) can bind one of its two heads (motor domains) to an empty lattice site vilfan2001a . The binding rate $`k_{+1}c`$ for this process is proportional to the solution concentration $`c`$ of the dimeric proteins. Successively, the dimer may either dissociate from the protofilament with a rate $`k_1`$ or also bind its second head to an unoccupied site in front (f) of or behind (b) the already bound head (with front we always refer to the direction pointing towards the “+” end of the microtubule, which is the direction of motion for most motors from the kinesin family). Since kinesin heads and microtubules are both asymmetric structures the corresponding binding rates $`k_{+2}^f`$ and $`k_{+2}^b`$ are in general different from each other. The reverse process of detaching a front or rear head occurs at rates $`k_2^f`$ and $`k_2^b`$. A reaction scheme with all possible processes and their corresponding rate constants is shown in Fig. 3.

In typical decoration experiments there is no external energy source, i.e. no ATP hydrolysis. Then the binding rates are not all independent of each other, but detailed balance dictates that the ratio of on- and off-rates has to equal the equilibrium binding constants

$$K_1=\frac{k_{+1}}{k_1},\text{and}K_2=\frac{k_{+2}^b}{k_2^b}=\frac{k_{+2}^f}{k_2^f}.$$
(1)
A particular coverage of the lattice is described as a sequence of dimers bound with both heads ($`D`$), one head only ($`1`$) and empty sites ($`0`$) (see Fig. 4). We denote the probabilities to find a certain lattice site in one of these states by $`2n_D`$, $`n_1`$ and $`n_0`$, respectively. Of course, normalization of the probabilities requires $`n_0+n_1+2n_D=1`$.

### II.1 Particle-hole symmetry

Important information about the steady state and the dynamics of the system can be gained already by exploiting the symmetries of the kinetic process. Fig. 5 reveals a “particle-hole” symmetry by showing a transformation that maps the system onto one with the same reaction scheme, albeit with transformed kinetic constants. The details of the transformation are listed in Table 1.

The “particle-hole” exchange operation replaces dimers bound on a single head (“particles”) by vacant sites (“holes”) and vice versa. Dimers bound with both heads are kept invariant under the transformation. To obtain a system with an equivalent reaction scheme the transition rates have to be transformed as well. This is achieved by exchanging the attachment and detachment rate of the first head, $`k_{+1}ck_1`$, and simultaneously switching the forward and backward binding rates of the second head, $`k_{2\pm }^fk_{2\pm }^b`$. As a consequence of this symmetry operations the equilibrium constant $`K_1c`$ is replaced by its reciprocal value while the binding constant of the second head $`K_2`$ is left unchanged.

This symmetry has immediate consequences on the equilibrium stoichiometry. Since the coverage in the steady state is solely a function of the quantities $`K_1c`$ and $`K_2`$ this symmetry implies that the mean number of dimers attached with both heads (which are mapped onto themselves) $`n_D(K_1c,K_2)`$ is invariant upon interchanging the attachment and detachment rates of the first head,

$$n_D(1/K_1c,K_2)=n_D(K_1c,K_2)$$
(2)
Similarly the mean total number of bound heads per lattice site (binding stoichiometry), $`\nu =2(n_1+n_D)`$, obeys the symmetry relation

$$\nu (1/K_1c,K_2)=2\nu (K_1c,K_2).$$
(3)
The symmetries can be best seen in a logarithmic-linear plot as shown in Fig. 6. From these relations we already conclude that the stoichiometry of dimers bound with both heads $`n_D`$ reaches its maximum at $`K_1c=1`$ and that the total stoichiometry is $`\nu =1`$ at that point.

### II.2 Experimental parameter values

The kinetic constants for the binding of kinesin on microtubules have been determined by several groups hackney94 ; ma95 ; moyer98 ; mandelkow98 . The binding constants for the first head in the presence of ATP and low ionic strength have the values moyer98  $`k_{+1}=20\mu \mathrm{M}^1\mathrm{s}^1`$ and $`k_1=25\mathrm{s}^1`$, leading to $`K_1=0.8\mu \mathrm{M}^1`$. $`k_1`$ is much smaller in the presence of AMP-PNP or in the absence of a nucleotide, about $`0.01\mathrm{s}^1`$ hancock99 . The binding constant $`K_2`$ can be estimated indirectly as the ratio between the detachment rate of the monomeric and the dimeric kinesin and has values between 2.7 (in the presence of ATP) and 20 (without a nucleotide) hancock99 . These values show that the model parameters depend strongly on the chemical conditions and the type of motor protein used in the experiment. A theoretical investigation allows for a systematic analysis and detailed classification of all the different regimes of binding kinetics in such a broad range of parameter values.

## III Equilibrium stoichiometry

In this section we review results obtained for the equilibrium stoichiometry of the dimer binding model vilfan2001a . In this section we will discuss the one-dimensional model where some short-range interaction between the dimers is taken into account; two-dimensional models are studied in Sec. V.

### III.1 Analytical solution for the binding stoichiometry

For the one-dimensional model the value of the mean occupation numbers in the steady state can quite easily be determined upon using detailed balance and the fact that the dimers have only a hard-core interaction. Because we have non-interacting dimers the probability to find a certain sequence of $`0`$’s, $`1`$’s and $`D`$’s (e.g. “0,1,D,D,1,D”, see Fig. 4) has to be invariant against permutations of these states. In such a random sequence the probability $`p_i`$ to find a particular state $`0`$, $`1`$ or $`D`$ at a certain site $`i`$ is given by

$$p_i=\frac{n_i}{n_0+n_1+n_D}.$$
(4)
Detailed balance requires that for each pair of possible configurations, their respective probabilities are in the same ratio as the transition rates between them. Hence the ratio of probabilities to find a sequence with a $`1`$ or $`0`$ at a certain place is

$$\frac{p_1}{p_0}=\frac{k_{+1}c}{k_1}=K_1c.$$
(5)
Similarly, we get for transitions between $`D`$ and $`01`$:

$$\frac{p_D}{p_0p_1}=\frac{k_{+2}^f}{k_2^f}=\frac{k_{+2}^b}{k_2^b}=K_2.$$
(6)
These two equations, together with the normalization condition, uniquely determine the values $`n_0`$, $`n_1`$ and $`n_D`$. The number of dimers bound with both heads is given by

$$n_D=\frac{1}{2}\frac{1}{2}\left(\frac{4K_1K_2c}{(1+K_1c)^2}+1\right)^{\frac{1}{2}}$$
(7)
and the number of dimers bound with one head

$$n_1=\frac{K_1c}{K_1c+1}\left(\frac{4K_1K_2c}{(1+K_1c)^2}+1\right)^{\frac{1}{2}}$$
(8)
The stoichiometry, i.e. the total number of heads per binding site $`\nu =2(n_1+n_D)`$ then reads

$$\nu =1+\frac{K_1c1}{K_1c+1}\left(\frac{4K_1K_2c}{(1+K_1c)^2}+1\right)^{\frac{1}{2}}.$$
(9)
The number of dimers bound with both heads per lattice site reaches its maximum $`n_D^{\mathrm{max}}=(11/\sqrt{K_2+1})/2`$ for $`K_1c=1`$. For an illustration of all these equations see Fig. 6.

### III.2 Alternative derivation of the binding stoichiometry

Here we present an alternative derivation of the binding stoichiometry which has the additional benefit that it can also be used for two-dimensional lattices where dimers can bind in either direction. The idea is to relate the equilibrium stoichiometry of the flexible dimer model to that of a stiff dimer model – a model where dimers can only bind as whole. Once we have specified the mapping one can determine the stoichiometry of the general model from the stoichiometry of the stiff dimer model in one dimension which has been known for quite a long time hill57  already and represents a special case of the solution by McGhee and von Hippel mcghee74

$$n=\frac{1}{2}\left(1\frac{1}{\sqrt{4K+1}}\right).$$
(10)
Note, that in this notation the binding constant $`K`$ contains the concentration $`c`$. In equilibrium, a mapping between the stiff and flexible dimer model is possible because of detailed balance and because transitions between states $`0`$ and $`1`$ are uncorrelated with the surrounding configuration. Let’s denote the equilibrium binding constant, the dimer and vacancy density in the stiff dimer model by $`K`$, $`\overline{n}_D`$ and $`\overline{n}_0`$, respectively. Then, detailed balance gives the following relation between these densities and the equilibrium binding constant, $`K=\overline{n}_D/\overline{n}_0^2`$. In the same way detailed balance implies $`K_2=n_D/(n_0n_1)`$ and $`K_1c=n_1/n_0`$ for the general dimer model. The absence of correlations allows us to group together the states $`1`$ and $`0`$ of the general model into a single state which we identify with the vacant state $`\overline{0}`$ of the stiff dimer model. This gives the mapping as summarized in Table 2.

From this table we can infer that the number of fully bound dimers per lattice site in the flexible dimer model is given as

$$n_D(K_1c,K_2)=\overline{n}_D(K)=\overline{n}_D\left(\frac{K_2K_1c}{(1+K_1c)^2}\right),$$
(11)
and the number of dimers bound with one head as

$$\begin{array}{c}n_1(K_1c,K_2)=\frac{K_1c}{1+K_1c}\overline{n}_0(K)=\hfill \\ \hfill =\frac{K_1c}{1+K_1c}\left(1\frac{1}{2}\overline{n}_D\left(\frac{K_2K_1c}{(1+K_1c)^2}\right)\right).\end{array}$$
(12)
In the one-dimensional case we can insert Eq. (10) and reproduce Eqs. (7) and (8). Note that the mapping reduces the number of parameters by one leading to a quite significant simplification. Though the stiff dimer model has not yet been solved analytically in two dimensions, the model with a single parameter can easily be solved by Monte-Carlo simulation (e.g., phares93b ). The entropy of a two-dimensional lattice fully covered with dimers is even known exactly kasteleyn61 ; temperley61 .

## IV Dynamics

In this section we will learn that the relaxation towards the equilibrium state described in the previous section is by no means a simple process characterized by a single time scale but actually shows quite interesting anomalous kinetics. We find different scenarios depending on the value of the equilibrium binding constants $`K_2`$ of the second head. If $`K_21`$ there are comparatively few double-bonded dimers in the steady state. Consequently cooperative effects and correlations are less important implying that the equilibration process is exponential with the time scale given by the rates $`k_{+1}c`$ and $`k_1`$. The dynamics changes drastically if the second head is very likely to bind, $`K_21`$. Then, the equilibrium state will consist mainly of dimers bound with both heads aligned with twice the period of the lattice. Their positions are correlated at a length scale given by the average distance between defects in the periodic pattern (either vacancies or single-bonded dimers). As we will detail in this section these correlations lead to a dynamics which is slowed down drastically as compared to the kinetics of individual dimers.

In the following we will mainly concentrate on the case $`K_1c1`$, to which we refer as the “stiff dimer model”, reflecting the fact that dimers are only likely to be found in the state with both heads bound. In this limit we can discuss the essential features of our model, which include the anomalous relaxation kinetics. Because of the particle-hole symmetry (Sect. II.1) which transforms $`K_1c1/(K_1c)`$, the results apply equally well for $`K_1c1`$. The other interesting limiting case is the situation with $`K_21`$ and $`K_1c1`$. In this limit the dynamics is essentially similar to the case $`K_1c1`$, but more complicated in detail because we have to consider all three states, i.e. empty sites, dimers bound with one head and dimers bound with both heads.

### IV.1 The stiff dimer limit

It turns out that most interesting features of the model can already be discussed if we restrict ourselves to the limit

$$K_21,K_1c1$$
(13)
with $`K=K_1K_2c`$ fixed at a value which allows the equilibrium state to consist mainly of dimers bound with both heads. The first condition of Eq. (13) states that the second head is unlikely to be found in the unbound state if it has a place where it can bind. The second condition implies that isolated single-site vacancies are unlikely to be occupied with single-bonded dimers. Therefore, our model reduces to a dimer deposition-evaporation model where the dimers can only bind and unbind with both heads at the same time. We will refer to this model as “stiff dimer model” in the following. The vacancy concentration $`n_0`$ in the steady state, given as $`1n_12n_D`$ in Eq. (7) and (8) then simplifies to

$$n_0=1/\sqrt{1+4K}.$$
(14)
The stiff dimer model describes dimers that bind at once to a pair of empty lattice sites and is characterized by the effective attachment rate of a whole dimer to an empty pair of sites, $`k_+c`$, an effective detachment rate $`k_{}`$, and a rate with which a dimer can move diffusively without detaching, $`r_d`$. Depending on whether the diffusion or the detachment rate is larger, we will distinguish between the *diffusive* and the *non-diffusive* case. The stiff dimer model is similar to conventional dimer deposition models privman92 , however with the difference that the attachment and detachment constants $`k_+`$ and $`k_{}`$ also depend on whether the neighboring sites are occupied or not (but we will show that this dependence becomes irrelevant in the non-diffusive case).

#### IV.1.1 Mapping between the general model and the stiff dimer model

In the following we describe how the kinetics of the stiff dimer model follows from our general dimer model. In the limit we are considering, the detachment of a single-bound dimer is always much faster than the attachment of a new one ($`k_1k_{+1}c`$, following from $`K_1c1`$) and the attachment of the second head much faster than its detachment ($`k_{+2}^{f,b}k_2^{f,b}`$, following from $`K_21`$). We do not put any other limitations on the reaction rates for now.

In a general treatment, the attachment and detachment rates also depend on the occupancy of the neighboring sites. We therefore have to distinguish between four cases: both neighbors empty, front neighbor occupied, rear neighbor occupied and both neighbors occupied.

* Case (i): Both neighbors empty. A dimer attaches with its first head at the rate $`k_{+1}c`$. Afterwards, the second head can attach at the total rate $`k_{+2}^f+k_{+2}^b`$, while the first head can detach with rate $`k_1`$. This leads to the total attachment rate per lattice site

$$k_+^Nc=\frac{k_{+1}c(k_{+2}^f+k_{+2}^b)}{k_1+k_{+2}^f+k_{+2}^b}.$$
(15)
The process of detachment starts with the detachment of the front or rear head. The detachment gets completed if the other head detaches too (rate $`k_1`$) and fails if the detached head reattaches again (rate $`k_{+2}^f+k_{+2}^b`$). The detachment rate then becomes

$$k_{}^N=\frac{k_1(k_2^f+k_2^b)}{k_1+k_{+2}^f+k_{+2}^b}.$$
(16)
* Case (ii): Front neighbor occupied. There are two pathways that lead to binding next to an occupied site. The first head can either attach on the site next to the occupied one or one site further away, both at the rate $`k_{+1}c`$. In the first case the probability that the second head will attach before the first one detaches is $`k_{+2}^b/(k_1+k_{+2}^b)`$. In the second case the probability that the second head attaches in front of the first one is $`k_{+2}^f/(k_1+k_{+2}^f+k_{+2}^b)`$. These two terms taken together give the attachment rate to the given pair of sites

$$k_+^Mc=k_{+1}c\left(\frac{k_{+2}^b}{k_1+k_{+2}^b}+\frac{k_{+2}^f}{k_1+k_{+2}^f+k_{+2}^b}\right).$$
(17)
If the process of detachment starts with detaching the front head (rate $`k_2^f`$), the probability that the whole dimer will detach afterwards is $`k_1/(k_1+k_{+2}^f+k_{+2}^b)`$. If it starts with detaching the rear head (rate $`k_2^b`$), this probability is $`k_1/(k_1+k_{+2}^b)`$. This gives the total detachment rate

$$k_{}^M=k_1\left(\frac{k_2^b}{k_1+k_{+2}^b}+\frac{k_2^f}{k_1+k_{+2}^f+k_{+2}^b}\right).$$
(18)
* Case (iii): Rear neighbor occupied. This case is analogous to the previous one, except that the indices $`f`$ and $`b`$ have to be interchanged. The attachment rate becomes

$$k_+^Pc=k_{+1}c\left(\frac{k_{+2}^b}{k_1+k_{+2}^f+k_{+2}^b}+\frac{k_{+2}^f}{k_1+k_{+2}^f}\right)$$
(19)
and the detachment rate

$$k_{}^P=k_1\left(\frac{k_2^b}{k_1+k_{+2}^f+k_{+2}^b}+\frac{k_2^f}{k_1+k_{+2}^f}\right).$$
(20)
* Case (iv): Both neighbors occupied. In this case the dimer binds into a gap of exactly two sites. If the first head attaches to the front site in the gap (rate $`k_{+1}c`$), the probability for the second one to attach before the first one detaches is $`k_{+2}^b/(k_1+k_{+2}^b)`$. Together with the analogous case in which the first head attaches to the rear site inside the vacancy, the total attachment rate becomes

$$k_+^Dc=k_{+1}c\left(\frac{k_{+2}^b}{k_1+k_{+2}^b}+\frac{k_{+2}^f}{k_1+k_{+2}^f}\right).$$
(21)
After the front head has detached (rate $`k_2^f`$), the probability that the whole dimer will follow is $`k_1/(k_1+k_{+2}^f)`$. Adding the pathway starting with the detachment of the rear head gives the detachment rate

$$k_{}^D=k_1\left(\frac{k_2^b}{k_1+k_{+2}^b}+\frac{k_2^f}{k_1+k_{+2}^f}\right).$$
(22)

Of course, all these rates obey detailed balance, which states that

$$\frac{k_+c}{k_{}}=K=K_1K_2c.$$
(23)
Processes in which one head detaches on one side and subsequently attaches on the other side also lead to an explicit diffusion of attached dimers. A diffusive step forwards occurs if the rear head detaches (rate $`k_2^b`$) and reattaches on the front side (probability $`k_{+2}^f/(k_1+k_{+2}^f+k_{+2}^b)`$). The hopping rate (which is, of course, equal for backward steps) then reads

$$r_d=\frac{k_{+2}^fk_{+2}^b}{K_2(k_1+k_{+2}^f+k_{+2}^b)}.$$
(24)
So far we have shown that in the limit $`K_21`$ our model can be mapped onto a simple dimer model. But, depending on parameters, this can happen in two qualitatively different ways. First, if the detachment rate is smaller than the diffusion rate $`k_{}r_d`$, which is equivalent to the condition

$$k_1\frac{1}{\frac{1}{k_{+2}^f}+\frac{1}{k_{+2}^b}},$$
(25)
we get the diffusive limit. In the opposite limit, a dimer is unlikely to make a diffusive step before it detaches and we call that the non-diffusive limit. In the non-diffusive limit, the effective dimer attachment ($`k_+c`$) and detachment ($`k_{}`$) rates become equal in all four cases, i.e. they become independent of the surroundings of a dimer.

To summarize, the stiff dimer model in the non-diffusive limit, to which we will concentrate in the following, describes dimers that can bind to any free pair of lattice sites with equal rate, namely $`k_+c`$ and can detach with a rate $`k_{}`$. The dimer cannot move along the lattice and it always stays joined (cannot regroup its two parts with parts of other dimers). Although dimers cannot diffuse directly, in the case of high lattice coverage we can still introduce an effective diffusion which describes dimers near a gap detaching and reattaching on the other side of the gap. This simplified model still captures all the essential aspects of the non-equilibrium dynamics.

#### IV.1.2 Related dimer models

Similar dimer adsorption models have been studied previously privman92 ; barma93 ; stinchcombe93 ; eisenberg97 . Privman and Nielaba privman92  studied the effect of diffusion on the dimer deposition process (neglecting dimer desorption). There are several key differences to our model, the most important of which is that diffusion without detachment results in a 100 % saturation coverage, whereas a model with detachment leads to a limiting coverage whose value depends on the binding constants of the first and second head (see Eq. 9). This also has important implication on the dynamics as discussed below. For example, as a consequence of a finite coverage the final approach to equilibrium is not a power law but exponential and there are, in addition, interesting temporal correlations in the fluctuations in the steady state.

Stinchcombe and coworkers barma93 ; stinchcombe93  studied the effect of detachment on the adsorption kinetics but allowed for regrouping of attached dimer molecules (two monomers that belonged to different dimers during attachment can form a dimer and detach together). While such processes are allowed for some types of inorganic dimers, they are certainly forbidden for dimer proteins like kinesin where the linkage between its two heads is virtually unbreakable. If regrouping is allowed the steady state auto-correlation functions for the dimer density shows an interesting power-law decay $`t^{1/2}`$ barma93 ; stinchcombe93 . This behavior can be directly linked to the gapless spectrum in an associated spin model stinchcombe93 . If regrouping is forbidden this power-law decay is lost (and becomes an exponential to leading order) due to the permanent linkage between the two heads of the dimer. Intuitively this may be understood as follows: only if regrouping of dimers is allowed are there locally jammed configurations (Néel-like states, “101010”, with alternating occupied and unoccupied sites in which neither attachment nor detachment of dimers is possible) in the final steady state which slow down the dynamics. The difference between the occupation numbers on even and odd sites then represents a conserved quantity slowing down the dynamics. In our model, the conservation law is trivial since the difference between occupation numbers on even and odd sites always disappears.

In addition, we will see in Sec. IV.1.4 that the autocorrelation functions of the dimer and the vacancy occupation number show strong differences in shape and typical times scales of relaxation.

#### IV.1.3 Numerical analysis of the non-diffusive stiff dimer model

To study the kinetics we choose the initial condition as typically used in an experiment, namely an empty lattice. Figure 8 shows simulation data for the average vacancy concentration as a function of time for a set of binding constants $`K=k_+c/k_{}K_1K_2c`$. We find qualitatively very different approaches to the final steady state depending on the value of the binding constant $`K`$. For $`K1`$, where the off-rate $`k_{}`$ is much larger than the on-rate $`k_+c`$, there is no crowding on the lattice and the dimeric nature of the molecules does not affect the approach to equilibrium. Hence one gets, like for monomers, an exponential behavior with a decay rate $`k_{}`$. In the opposite limit, $`K1`$, there is a pronounced two-stage relaxation towards the steady state. The vacancy concentration as a function of time reveals four regimes, an initial attachment phase, followed by an intermediate plateau, then a power-law decay and finally an exponential approach towards equilibrium.

##### Initial binding

At short time scales, $`tk_{}^1`$ and $`tr_d^1`$, when only deposition processes are frequent but detachment or diffusion processes are still very unlikely, the kinetics of the model belongs to the class of problems frequently referred to as random sequential adsorption (RSA) models (reviewed by Evans evans93  and Talbot *et al.* talbot00 ). RSA Models have among others been used to describe the chemisorption of inorganic dimers like $`\mathrm{O}_2`$ on crystal surfaces, the binding of reagents to organic polymers or chemical reactions on the polymer chain.

The deposition of stiff dimers has first been studied by Flory flory39 . Later, it was pointed out by Page page59  (see also a note by Downton downton61 ) that the final distribution depends on whether dimers bind with both ends in parallel (“*standard model*”) or with one end first (“*end-on model*”). A more recent article addressing this point and giving analytical solutions for both cases was published by Nord and Evans nord90 . In the “standard model” the dimers bind with equal rates to all free pairs of lattice sites. In the “end-on” the first end of a dimer binds with equal rates to any free site. The second end then binds to either side if there is space, or to the only available side if the other one is occupied. This has the consequence that the probability that a dimer binds to either end of a gap is $`1.5`$ times higher than that it binds to a certain pair of sites in the middle of the gap. This obviously leads to some clustering of dimers reduces the number of gaps that remain in the end.

It is even possible to exactly calculate the time-dependence of the lattice coverage mcquistan68 . In the “standard model” the vacancy density during the RSA phase follows mcquistan68 ; evans93

$$n_0(t)=\mathrm{exp}\left(2+2e^{k_+ct}\right).$$
(26)
The vacancy concentration locks at an intermediate plateau

$$n_0^{\mathrm{plateau}}=e^20.1353$$
(27)
This Flory plateau represents a configuration in which all remaining vacancies are isolated, causing the system to be unable to accommodate for the deposition of additional dimers. In the “end-on” model the plateau vacancy concentration is

$$n_0^{\mathrm{plateau}}=\sqrt{2\pi e}(\mathrm{Erf}(\sqrt{2})\mathrm{Erf}(1/\sqrt{2}))10.1233,$$
(28)
slightly smaller than in the non-diffusive limit.

When applying the RSA results to our model, we have to distinguish between the diffusive and the non-diffusive case. In the *non-diffusive* case, the binding rate at all pairs of sites, regardless whether their neighboring sites are empty (15), one of them (17, 19) or both (21) are occupied, is equal. This corresponds to the “standard model”. The time-dependent gap density is given by (26) and its plateau value by (27).

In the *diffusive* case, the binding rates on both ends of an interval, (17) and (19), add up to a rate three times as large as the binding rate in the middle of an interval (15), while the attachment rate for a pair of sites with both neighbors occupied (21) is twice as large. Therefore, the RSA phase in the diffusive limit corresponds to that of the “end-on model”. The plateau gap density is then given by (28). In both cases, the final configuration is independent of the asymmetry in the binding rates.

##### Power-law regime

The secondary relaxation process towards the final equilibrium state is much slower than the initial random sequential adsorption process. Starting out of a jammed configuration in the Flory plateau the dynamics shows a broad time domain with a power-law $`t^{1/2}`$ instead of a simple exponential decay. Similar multi-stage relaxation processes have been observed in dimer models with diffusion but no detachment privman92 , the key difference being that the detachment process implies that the steady state has a finite vacancy density and the final approach to the steady state remains not a power law but becomes exponential. There are also interesting similarities. In particular, in both models a large portion of the final approach to the steady state is mediated by the annihilation of vacancies. This behavior can be explained by introducing a particle representation in the following way (analog to the adsorption-diffusion model privman92 ). We denote each vacancy on the lattice as a “particle” $`A`$, and each bound dimer as an inert state ($`00`$). The detachment of a dimer then corresponds to a pair creation process $`00AA`$, and the decoration process to pair annihilation $`AA00`$. Since we consider the limit $`K1`$, states with two vacancies $`A`$ on neighboring sites have a very short lifetime. We may thus introduce a coarse-grained model by eliminating these states (see Fig. 9).

Then processes like $`A00AAA00A`$ result in an effective diffusion for particle $`A`$ with a hopping rate $`r_{\mathrm{hop}}=k_{}/2`$ ($`k_{}`$ is the rate of the first transition and $`1/2`$ the branching probability for the second) and an effective step width of two lattice sites. The hopping rate is different for the last step before two particles annihilate ($`A00A00AA`$), namely $`k_{}/3`$. Pair annihilation, $`AA00`$, occurs with a rate $`k_+c`$, or $`2k_{}/3`$ from the state $`A00A`$. Pair creation, on the other hand, occurs mainly through the process $`000000AAAAAAA00A`$; the corresponding rate is $`k_{}/3K`$ per lattice site and hence largely suppressed with respect to the annihilation process as long as the particle concentration is far above its steady-state value. These processes are summarized in Tab. 3. Other processes involve more particles and are of higher order in terms of a power series in $`K^1`$. They are negligible if the number of vacancies is small, $`n_A1`$, and the binding constant is large, $`K1`$. In summary, for $`K1`$, our dimer model can be mapped onto a one-particle reaction-diffusion model $`A+A0`$. Pair creation processes, $`0A+A`$, are highly suppressed during the first stage of the relaxation process. The give a significant contribution to the dynamics close to the steady state where the relaxation becomes exponential instead of algebraic.

Models like the simple reaction-diffusion model $`A+A0`$ show interesting non-equilibrium dynamics privman97 ; mattis98 ; schuetz98 . One can show torney83 ; toussaint83  that asymptotically the particle density $`n_A(t)`$ decays algebraically, $`n_A(t)t^{1/2}`$, which nicely explains the slow decay observed in simulation data (see Fig. 8). Note that a mean-field like rate equation

$$_tn_A(t)n_A^2(t)$$
(29)
would predict $`n_A(t)t^1`$. In our analysis we can even go beyond the asymptotic scaling analysis and try to compare with exact solutions of the model for a random initial distribution with density $`p`$ by Krebs et al. krebs95 . They find (adapted to our situation with two-site hopping)

$`n_A(t)={\displaystyle \frac{1}{2\pi }}{\displaystyle \underset{0}{\overset{2}{}}}𝑑u{\displaystyle \frac{\sqrt{u(2u)}\left(n_0^{\mathrm{plateau}}\right)^2e^{16ur_{\mathrm{hop}}t}}{u\left(u(\frac{1}{2}n_0^{\mathrm{plateau}})+\left(n_0^{\mathrm{plateau}}\right)^2\right)}}.`$ (30)
Its asymptotic limit (first determined by Torney and McConnel torney83 ) reads

$$n_A=(32\pi r_{\mathrm{hop}}t)^{1/2}$$
(31)
Note that it is independent of the initial particle concentration in the Flory plateau $`n_0^{\mathrm{plateau}}`$.

Our Monte-Carlo data (see Fig. 8) are in excellent agreement with the predictions of Eq. 30. Minor deviations at times between the plateau and the power-law decay are due to the assumption of a random particle distribution underlying the derivation of Eq. 30. In reality the state after initial binding (RSA) contains some correlations which however do not reach beyond a vacancy (Markov-shielding) evans93 . In other words, the sizes of occupied areas are not exponentially distributed, but they are uncorrelated to each other and the higher-order correlation functions decay as fast as 2-point correlations (super-exponentially). Short ranged correlations in the initial configuration do have some effect on concentration at intermediate times and this causes the deviation between the simulation data and the theoretical curve which, however, is not large. According to family91 ; santos97  correlations can affect the concentration in the asymptotic limit even if the pair correlation function is short ranged and the order-$`n`$ correlation functions decay at a length-scale of $`𝒪(n)`$, but this is not the case in our model.

The results from the $`A+A0`$ model also become invalid for very long times where the particle concentration comes close to its equilibrium value. In this limit the dynamics becomes scale-invariant racz85 . The particle concentration can be written in scaling form

$$n_A(t)=K^{1/2}\widehat{n}(k_{}t/K)$$
(32)
with a diverging characteristic time scale $`\tau _KK`$.

#### IV.1.4 Steady-state autocorrelation functions

A limitation of the above mapping becomes evident if one considers the equilibrium autocorrelation functions. Contrary to conventional $`A+A0`$ models there are three different autocorrelation functions with different functional forms (Fig. 10). In the following we will use $`\widehat{n}_{(i,i+1)}`$ as the occupation number (0 or 1) of a dimer on the pair of sites $`(i,i+1)`$ and $`\widehat{n}_{0_i}`$ as the occupation number of a vacancy at the site $`i`$ (1 if there is a vacancy, 0 otherwise). The dimer autocorrelation function

$$C_D(\tau )=\widehat{n}_{(i,i+1)}(t)\widehat{n}_{(i,i+1)}(t+\tau )\widehat{n}_{(i,i+1)}(t)^2$$
(33)
describes the probability to find a dimer at a given pair of sites simultaneously at times $`t`$ and $`t+\tau `$. Similarly we define the vacancy autocorrelation function

$$C_0(\tau )=\widehat{n}_{0_i}(t)\widehat{n}_{0_i}(t+\tau )\widehat{n}_{0_i}(t)^2$$
(34)
describing the probability that the site $`i`$ is vacant at times $`t`$ and $`t+\tau `$. To calculate their values at $`\tau =0`$, we need the expectation values

$$\widehat{n}_{(i,i+1)}(t)=\widehat{n}_{(i,i+1)}^2(t)=n_D=\frac{1}{2}\left(1\frac{1}{\sqrt{4K+1}}\right)$$
(35)
and

$$\widehat{n}_{0_i}(t)=\widehat{n}_{0_i}^2(t)=n_0=\frac{1}{\sqrt{4K+1}}.$$
(36)
A third way of defining an autocorrelation function is to ask for the probability that a vacancy is either on the site $`i`$ or on the site $`i+1`$

$$\begin{array}{c}C_P(\tau )=<(\widehat{n}_{0_i}+\widehat{n}_{0_{i+1}}\widehat{n}_{0_i}\widehat{n}_{0_{i+1}})(t)\times \hfill \\ \hfill \times (\widehat{n}_{0_i}+\widehat{n}_{0_{i+1}}\widehat{n}_{0_i}\widehat{n}_{0_{i+1}})(t+\tau )>\\ \hfill (\widehat{n}_{0_i}+\widehat{n}_{0_{i+1}}\widehat{n}_{0_i}\widehat{n}_{0_{i+1}})(t)^2,\end{array}$$
(37)
with

$$\widehat{n}_{0_i}+\widehat{n}_{0_{i+1}}\widehat{n}_{0_i}\widehat{n}_{0_{i+1}}=\frac{2n_0}{1+n_0}=\frac{2}{\sqrt{4K+1}+1}.$$
(38)
The important difference between $`C_P`$ and $`C_0`$ originates from the fact that the vacancies diffuse by hopping over two lattice sites. This is why a vacancy diffusing along odd sites will not influence the function $`C_0`$ if the original vacancy was on an even site. The function $`C_P`$ does not make this distinction.

For $`K1`$ the autocorrelation functions become scale invariant as well. Their scaling form reads $`\widehat{C}_0(\widehat{t})=KC_0(t)`$, $`\widehat{C}_D(\widehat{t})=C_D(t)`$ and $`\widehat{C}_P(\widehat{t})=KC_P(t)`$ with $`\widehat{t}=k_{}t/K`$. The latter corresponds to the autocorrelation function in a reaction-diffusion model $`A+A0`$, which has recently been calculated analytically by Bares and Mobilia bares99

$$\begin{array}{c}\widehat{C}_P(\widehat{t})=\left(\frac{e^{2\widehat{t}}}{\sqrt{2\pi \widehat{t}}}\text{Erfc}\sqrt{2\widehat{t}}\right)\text{Erfc}\sqrt{2\widehat{t}},\hfill \\ \hfill \mathrm{with}\text{Erfc}(x)=\frac{2}{\sqrt{\pi }}_x^{\mathrm{}}e^{y^2}𝑑y.\end{array}$$
(39)
Note that the described solution is derived for a parameter set which requires a special relation between diffusion, pair-creation and annihilation rate to be fulfilled. In our model this is not the case. But, this difference becomes irrelevant in the scaling limit, since their and our model can be mapped onto each other by introducing a short-ranged interaction between particles. See also Ref. park01  for a comment about the validity of the analytical approximation used in Ref. bares99 . The observed deviations however do not affect the scaling limit.

The other two functions decay on the same time-scale, but with different prefactors. The reason is that even if a pair of vacancies annihilates, the system still keeps memory on whether the surrounding dimers were located on even or odd locations and this gives those correlation functions that distinguish between even and odd sites a longer decay time.

### IV.2 Two-particle model

Another interesting case arises if we consider the limit $`K_2\mathrm{}`$ with $`K_1c`$ being of the order of magnitude $`1`$. In this case neither the dimers bound with one head nor the empty lattice sites dominate and we have to introduce two particle types. As we are used to, particles $`A`$ should denote vacancies on the lattice. In addition, we introduce particles $`B`$, representing dimers bound with only one head.

Reactions $`AB`$ and $`BA`$ occur at rates $`k_{1+}c`$ and $`k_1`$. If we assume that both attachment rates of the second head $`k_{2+}^f`$ and $`k_{2+}^b`$ are significantly higher than $`k_{1+}c`$ and $`k_1`$, we are again dealing with the *diffusive model* (25). The diffusion mechanism for vacancies (particles A) is exactly the same as in the model with stiff dimers and according to (24) the hopping rate is given by Eq. (24). The diffusion rate of single-bonded dimers (particles B) has to be equal. This is due to the symmetry of our model upon exchanging particles A and B and the transition rates $`k_{1+}c`$ and $`k_1`$ (Tab. 1) and because the latter are irrelevant for the hopping rate. If a particle A and a particle B reach neighboring sites, they annihilate quickly, while two particles of the same type do not interact. To summarize, we obtain a reaction diffusion model of the type:

$`A`$ $`B`$
$`A+B`$ $`0`$
While its dynamics is more complex at short times, the model becomes equivalent to the $`A+A0`$ if the transitions between $`A`$ and $`B`$ become faster than the typical annihilation time, which is given as the diffusion time between two particles, displaced by the average distance between particles on the lattice, $`\frac{1}{r_{\mathrm{hop}}n^2}`$. Therefore the power-law behavior, given by (31) remains untouched by the fact that we are dealing with two different particle types. An example of a simulation in this regime is shown in Fig. 11.

### IV.3 Dynamics of interacting dimers in one dimension

In many biological systems the interaction between attached molecules plays an important role. For example, in the case of kinesin, there has been an observation of coexisting empty and decorated domains which can only be explained by an attractive interaction vilfan2001a . Similar observation has been made on actin decorated with double-headed myosin orlova97 . An interaction between the dimers can be introduced by assuming that a dimer is more likely to bind to a pair sites if one or two neighbors are already bound. The binding rate then becomes $`k_+^P=k_+^N=Ak_+^N`$ (one neighbor bound) or $`k_+^D=A^2k_+^N`$ (both neighbors bound). Similar, we assume that a dimer with one bound neighbor dissociates with rate $`k_{}^P=k_{}^M=Bk_{}^N`$ and that a dimer with two bound neighbors dissociates with rate $`k_{}^D=B^2k_{}^N`$. The equilibrium state of this model is still exactly solvable mcghee74 , but the interaction changes both relaxation stages quantitatively. First, the vacancy concentration on the intermediate plateau lowers since the interaction improves the formation of contiguous clusters during the first stage. The RSA phase can then be described in terms of Kolmogorov’s grain-growth model kolmogorov37 . The vacancy concentration after the initial binding is given as $`n_0\frac{1}{2}\sqrt{\frac{\pi r_{\mathrm{nuc}}}{2r_{\mathrm{growth}}}}`$ gonzalez74 ; evans83 , where $`r_{\mathrm{nuc}}`$ is the nucleation rate per free lattice site (in the simplest case, when $`K1`$, it is simply $`k_+`$) and $`r_{\mathrm{growth}}`$ the sum of speed with which both boundaries of a nucleus spread over the lattice (normally $`k_+^P+k_+^M=2Ak_+^N`$). The second effect of the interaction is that the diffusional relaxation slows down by the factor $`B`$ since the detachment rate decreases. And finally, the equilibrium vacancy concentration decreases mcghee74 . Nevertheless, interacting models show the same two-stage relaxation behavior. An example of a model with interaction is shown by the dot-dashed line in Fig. 12.

## V Dimer model with interaction

As we mentioned in Section IV.3, interactions between bound dimers can often play an important role. In the kinesin-microtubule system, evidence for the existence of an attractive interaction comes from two key experimental observations. The first one is that bound kinesin dimers are found to form two-dimensional crystalline lattices thormaehlen98 . Within the non-interacting model one can explain longitudinal correlations along a single protofilament, but by no means the alignment of dimers on different protofilaments. The second observation supporting an attractive interaction between bound dimers is the coexistence of empty and decorated microtubules in one and the same experiment vilfan2001a . Such a phase segregation is a clear indication that the system is in a coexistence regime with the strength of the dimer interaction above some critical value. The detailed nature of the interaction between the dimers is not yet known. It could be a direct interaction between kinesin heads and tails or some indirect interaction mediated through distortions of the underlying tubulin lattice or a combination of all of these possibilities. All of these mechanisms are consistent with the observation that the interaction is stronger on flat tubulin sheets than on cylindrical microtubules vilfan2001a .

Even if one assumes that there is only nearest-neighbor interaction, a general description would still require $`12`$ different coupling strengths ($`2`$, $`4`$ and $`6`$ for describing the attraction between two single-bonded dimers, two double-bonded dimers and between a single- and a double-bonded dimer, respectively). Upon neglecting that the tubulin lattice is actually rhombic and not orthogonal one can reduce the number of different coupling constants to $`9`$. These are still too many for a general discussion. Hence, mainly for simplicity, we will restrict ourselves to a model with only two interaction constants, $`J_L`$ acting in longitudinal and $`J_T`$ acting in transverse direction, as illustrated in Fig. 13. We assume that the interaction between two double-bonded dimers on neighboring protofilaments is zero if they are displaced by one lattice site or more (Fig. 13d). Otherwise, the transverse interaction is assumed to have the strength $`J_T`$ between two double-bonded dimers (Fig. 13c) and $`\frac{1}{2}J_T`$ between a double- and a single-bonded dimer (Fig. 13b). We use the same notation as in Sec. IV.3 and introduce factors $`A_i`$ and $`B_i`$ measuring the change in the attachment and detachment rate by the presence of a neighboring dimer (with $`i`$ indicating the relative position of the neighbors as shown in Fig. 13),

$`k_{1,2+}^{\text{with neighbour}}`$ $`=A_ik_{1,2+}^{\text{without neighbour}},`$
$`k_{1,2}^{\text{with neighbour}}`$ $`=B_ik_{1,2}^{\text{without neighbour}}.`$
Detailed balance states that

$$J_i=k_BT\mathrm{ln}\frac{A_i}{B_i}.$$
(40)
For simplicity we further assume that the binding and the unbinding rate are always changed by the same factor, $`A_i=1/B_i`$; this assumption is irrelevant in equilibrium, but it simplifies the kinetics.

Fig. 14 shows the stoichiometry curves, equivalent to those in Fig. 6, for the interacting model. The most dramatic effect of the interaction is that it broadens the plateau where most dimers are double-bonded, giving a stoichiometry of one head per lattice site. If the interaction is strong enough (see e.g. the curve with $`J/k_BT=2`$ in Fig. 14), the stoichiometry curves show quite a steep rise from zero to the plateau at $`\nu =1`$. This is quite indicative of a phase transition, albeit smeared somewhat out due to finite-size effects.

To determine the critical point we make two further simplifications. First, we assume that we are in the stiff dimer limit, $`K_21`$, as we did in Sec. IV.1. This stiff dimer model can be interpreted as a spin model with *axial next nearest neighbor interaction*, also known as ANNNI models, reviewed by Selke selke88 . In terms of these ANNNI models, our model is a rather peculiar limit since it has an infinitely repulsive interaction between nearest neighbors and attractive interaction between next-nearest neighbors in longitudinal direction, as well as an attractive interaction between nearest neighbors in transverse direction. In the literature on ANNNI models mainly the opposite situation was studied since it leads to frustrations, while the “zero-temperature” state of our model is trivial, namely dimers aligned, let us say, the even sites on the protofilaments. This justifies an approximation which further simplifies the model, namely by assuming that dimers can only bind to even sites. Then the problem simplifies to a lattice-gas model with an effective Hamiltonian

$$=\underset{i,j}{}\left(J_L\widehat{n}_{2i,j}\widehat{n}_{2i+2,j}J_T\widehat{n}_{2i,j}\widehat{n}_{2i,j+1}\mu \widehat{n}_{2i,j}\right),$$
(41)
where the coordinates $`(2i,j)`$ run over all even lattice sites, $`\widehat{n}_{2i,j}`$ are the occupancy numbers ($`0`$ or $`1`$) on those sites and $`\mu =k_BT\mathrm{ln}(K_1K_2c)`$ is the chemical potential of adsorbed dimers. The lattice gas model can be mapped onto the 2D Ising model yeomans92  with Hamiltonian

$$^{}=\underset{i,j}{}\left(J_L^{}s_{i,j}s_{i+1,j}J_T^{}s_{i,j}s_{i,j+1}ϵs_{i,j}\right)$$
(42)
with exchange constants

$$J_T^{}=\frac{1}{4}J_T,J_L^{}=\frac{1}{4}J_L,$$
(43)
and spin variables

$$s_{i,j}=2\widehat{n}_{2i,j}1.$$
(44)
The spin $`s`$ has the value $`+1`$ if it is parallel and $`1`$ if it is antiparallel to the magnetic field $`ϵ`$,

$$ϵ=\frac{\mu +J_L+J_T}{2}.$$
(45)
The critical temperature of the two dimensional Ising model was first determined by Kramers and Wannier kramers41 . In the absence of an external field, the model has also been exactly solved by Onsager onsager44 . The condition for the critical point reads

$$\mathrm{sinh}\frac{J_L}{2k_BT}\times \mathrm{sinh}\frac{J_T}{2k_BT}=1,$$
(46)
or $`J_c=1.76k_BT`$ in the isotropic case ($`J_L=J_T`$).

The exact solution helps us to determine the critical point, i.e. the interaction strength needed to obtain a phase transition between an empty and a decorated state. The full form of the binding stoichiometry as a function of concentration, however, would be equivalent to calculating the magnetization of the two-dimensional Ising model in the presence of an external field, which has not yet been achieved analytically. Therefore, one has to rely on computer simulations to obtain the stoichiometry curves. An overview on Monte-Carlo simulations on lattice gas models can be found in Ref. binder89 . The situation simplifies when one is far away from the critical point, i.e. when the coupling is much stronger than its critical value, $`JJ_c`$. This corresponds to the low-temperature limit in the Ising analogue, where the spin polarization is saturated and points into the direction of the external field. The number of spins with orientation $``$ is given as $`n=N\mathrm{\Theta }(ϵ)`$, where $`\mathrm{\Theta }`$ is the Heaviside function. Back in the adsorption model, according to Eqs. (43)-(45), the condition reads

$$\nu =\mathrm{\Theta }(\mu +J_L+J_T).$$
(47)
This expression can also be understood directly. The boundary of a decorated domain will normally consist of many edge elements (each one with three out of four neighboring sites occupied) and some corner elements (with two out of four neighboring sites occupied). If the corner elements are stable, the decorated domain will grow; if they are unstable, it will shrink. The stability condition of a corner element with one bound longitudinal neighbor and one bound transverse neighbor is given by Eq. (47). Deviations from (47) are possible due to finite size effects. This is because small patches have a higher curvature in the boundary and are therefore less stable.

### V.1 Dynamics of the interacting 2D model

We have shown in Sec. IV.1.3 that the annealing of gaps between bound dimers can be slowed down substantially as compared to the transition rates for single dimers. This effect is even stronger in a two dimensional system. Here linear domain boundaries emerge instead of point defects (gaps). These boundaries can annihilate when two of them join. But their diffusion is significantly slower than the diffusion of single gaps, roughly inverse proportional to the number of sites in the boundary.

The nonequilibrium dynamics of decorating an initially empty two-dimensional lattice is an interesting and complex problem. It has first been studied by Becker and Döring becker35  (reviewed in Ref. lewis78 ). A crucial concept is the critical nucleus, a certain number of adsorbed atoms (here dimers) necessary to form a stable nucleus able to grow and spread over the whole lattice. An analytical expression for the nucleation time can be given if the interaction is strong enough, and the critical cluster size becomes as small as four, three or two molecules. For larger critical cluster sizes a rough estimate says that the nucleation rate is proportional to the kinesin concentration to a power that is typically half the number of nucleus-forming units.

There are different scenarios depending on the relative magnitude of the time to form a critical nucleus, $`t_{\mathrm{nuc}}`$, and its growth time, $`t_{\mathrm{growth}}`$, i.e. the time it takes such a nucleus to grow and cover the whole surface (see Fig. 15a).222For a Java applet showing the nucleation in the two-dimensional dimer model see http://www.ph.tum.de/~avilfan/decor/ . If the nucleation time is larger than the growth time,

$$t_{\mathrm{nuc}}t_{\mathrm{growth}},$$
(48)
the whole lattice will most probably be completely covered with dimers before other nucleation sites can form anywhere else on the lattice. The final state is a defect free lattice generated from a single nucleus. This is consistent with experimental observations thormaehlen98 . At the same time single-site nucleation also explains how empty and fully decorated microtubules can be found coexisting in one and the same experiment vilfan2001a . When the decoration starts, nuclei will form on certain microtubules and totally decorate them, until the concentration of free kinesin in solution drops below the critical value, where the equilibrium phase is the empty lattice. A numerical study of the nucleation process vilfan2001a  gives us an estimate of the minimum interaction strength needed to fulfill the condition (48) as $`J3k_BT`$. This strongly indicates that there is a interaction between kinesin dimers which is of the order of $`3k_BT`$.

If the nucleation time is smaller than the growth time,

$$t_{\mathrm{growth}}t_{\mathrm{nuc}},$$
(49)
there is a certain probability that a second critical nucleus forms on the lattice before the first one created had a chance to cover the whole lattice. The probability of such or higher order events grows with decreasing nucleation time. Then, the lattice will be covered by several domains, locally ordered areas which may be out of phase on the microtubule lattice. In order to finally reach a homogeneous decoration with only one predominant domain a secondary process is needed which leads to a coarsening of the initial grain-like structure. This secondary process includes domain wall wandering and annihilation induced by the detachment and re-attachment of kinesin dimers. As mentioned above, such annealing processes are expected to be much slower than the unhindered growth of individual domains. Fig. 15b summarizes the possible phases and dynamic regimes for interacting dimers decorating a two-dimensional lattice.

## VI Interacting $`k`$-mers in one dimension

Non-interacting $`k`$-mers ($`k>3`$) are well described with the mean-field kinetics of the model $`kA0`$ which predicts a gap density decay $`t^{1/(k1)}`$ nielaba92 ; privman92a . In the continuous limit, also called continuous sequential adsorption (CSA) or “car parking problem” the length of unoccupied space decays as $`1/\mathrm{log}t`$ in the mean-field regime krapivsky94 ; jin94 .

The dynamics becomes more intriguing when a strong attraction between the $`k`$-mers is introduced. An example of such a system is tropomyosin vilfan2001c , a protein that plays a key role in the regulation of muscle contraction. In muscle each tropomyosin molecule binding on an actin filament occupies seven actin monomers (although other isoforms covering 5 or 6 monomers exist as well hitchcock-degregori94 ) and thereby obstructs the action of myosin motors. The binding of tropomyosin to actin is strongly cooperative – the binding constant next to a bound molecule is 1000 times larger than at an isolated place wegner79 .

The initial binding phase is analog to the problem of interacting dimers and can again be described using a grain growth model kolmogorov37 ; evans93 ; vilfan2001c . After the initial phase, gaps with sizes between $`1`$ and $`k1`$ sites will remain. Their sizes are distributed randomly with equal probabilities of $`1/(k1)`$ for each gap size. What follows is again a process on a much longer time-scale in which $`k`$-mers at edges of the gaps can detach and re-attach. If they detach from one side of the gap and then re-attach at the other side, the gap moves by $`k`$ sites in one direction (see Fig. 16a).

The rate of such diffusive steps is given as the detachment rate of a molecule on e.g. the minus side of a block multiplied by the probability that a molecule re-attaches on the plus side of the neighboring block before another $`k`$-mer can re-attach to the position where the first one has detached:

$$r_{\mathrm{hop}}=k_{}^M\frac{k_+^P}{k_+^P+k_+^M}=k_{}^M\frac{k_{}^P}{k_{}^P+k_{}^M}=\left(\frac{1}{k_{}^P}+\frac{1}{k_{}^M}\right)^1.$$
(50)
Two neighboring gaps (see Figs. 16b and c), each containing between $`1`$ and $`k1`$ sites with randomly distributed probabilities, can either join to a single gap (Fig. 16b) or, if their total size exactly fits one $`k`$-mer, annihilate (Fig. 16c). If the original gap sizes are $`g_1`$ and $`g_2`$, the joined gap has size $`(g_1+g_2)modk`$. The annihilation will take place in $`k1`$ out of $`(k1)^2`$ cases, therefore its probability is $`1/(k1)`$. Any other gap size will be the outcome in $`k2`$ cases, therefore having the probability $`(k2)/(k1)^2`$.

Again we can represent gaps as particles $`A`$ hopping randomly along the lattice and joining or annihilating when two of them meet. The pair creation processes become relevant only if the particle concentration approaches its equilibrium value and will be neglected for now. Then, we can map our model to a diffusion-coagulation-annihilation model consisting of following reactions

$`A+A`$ $`A`$ with probability $`(k2)/(k1),`$
$`A+A`$ $`0`$ with probability $`\mathrm{\hspace{0.33em}1}/(k1).`$
The fundamental difference between the interacting $`k`$-mer and the dimer model is that the $`k`$-mer model not only contains pair annihilation but also pair coagulation processes, $`A+AA`$. According to Ref. lee94  all diffusion-annihilation models of the type $`2AlA`$ (or generally $`mAlA`$) belong to the same universality class. In the asymptotic limit, the “particle” concentration can then be borrowed from the exactly solvable $`A+A0`$ model lushnikov87 ; torney83 . Adapted to our $`k`$-mer model it reads

$$n_A(t)=\frac{2}{2l}\frac{1}{\sqrt{8\pi \overline{r}_{\text{hop}}t}}$$
(51)
with $`l=(k2)/(k1)`$ and $`\overline{r}_{\text{hop}}=k^2r_{\text{hop}}`$; the second relation results from the fact that in each diffusive step a gap jumps over $`k`$ sites. Hence, we finally obtain

$$n_A(t)=\frac{k1}{k^2\sqrt{2\pi r_{\text{hop}}t}}.$$
(52)
Interestingly, the asymptotic particle concentration is independent of its initial value, i.e. the final gap concentration at long times does not depend on the intermediate gap concentration $`n_A^0`$. Note that it is even independent of the solution concentration $`c`$.

Of course, the mapping to the $`A+AlA`$ model is valid only as long as the concentration of particles $`A`$ is well above its equilibrium concentration. Otherwise, particle splitting events like $`AA+A`$ also become relevant. If the concentration $`c`$ is high enough, the $`t^{1/2}`$ regime will be followed by a regime where, due to these splitting events, the average gap size will finally approach one site, which then are randomly distributed over the lattice. The number of gaps will decay according to the mean-field prediction $`t^{1/(k1)}`$, as it does in the non-interacting $`k`$-mer model nielaba92 ; privman92a . An example of this behavior is shown in Fig. 17. When the system approaches equilibrium, processes of pair creation $`0A+A`$ become important as well. Similar to the dimer model, these processes leads to a final exponential relaxation towards equilibrium. The equilibrium gap size and distribution of the interacting $`k`$-mer model are known exactly from a study by McGhee and von Hippel mcghee74 .

## VII Summary and Outlook

In this contribution we have discussed the kinetics of some macromolecular assembly processes relevant for the formation of functional structures in cells. In particular, we were interested in the dynamics of ligand-substrate binding, where the substrate is a one- or two-dimensional lattice and the ligands are dimers or oligomers. As examples we picked binding of dimeric kinesin on microtubules and tropomyosin on actin filaments. There are other related systems where protein ligands bind to biological macromolecules, e.g. proteins to DNA or antibodies to viruses.

Our theory allows us to study the effects of steric constraints, binding rates and interaction between neighboring proteins on the binding dynamics and binding stoichiometry. Our key results are as follows. The collective dynamics turns out to be not only much slower than the kinetics of single molecules but it also shows some interesting dynamic anomalies. Quite generally, we find that the binding kinetics goes through several qualitatively different stages. At small times when deposition processes are frequent but detachment processes are still unlikely, the kinetics can be described in terms of sequential deposition models. Such models have been analyzed in various contexts years ago so that there are many exact results available. After this initial phase the system is left in a “metastable” or “locked” configuration, i.e. a configuration with many small vacant regions between blocks of bound molecules prohibiting binding of additional ligands to the substrate. As a consequence one obtains an intermediate plateau regime for the binding stoichiometry. For escaping from this locked configuration vacant regions have to merge in order to accommodate for the deposition of additional ligands. These are extremely slow processes mediated by a sequence of detachment and re-attachment processes. Upon identifying the gaps with particles $`A`$ we have shown that the dynamics in this regime can be explained by mapping it onto reaction-diffusion models. For stiff dimers the corresponding model is $`A+A0`$ resulting in a gap density that decays with a power law $`t^{1/2}`$ before it approaches the final equilibrium value. Taking into account that actual biological dimers have some level of flexibility the model had to be generalized to include a second species of particles $`B`$, representing dimers bound with only one head. The corresponding reaction-diffusion model is then of the type $`AB`$ and $`A+B0`$. While its dynamics is more complex at short times, the model becomes equivalent to the $`A+A0`$ if the transitions between $`A`$ and $`B`$ become faster than the typical annihilation time. If, on the other hand, one introduces a strong diffusion of adsorbed dimers the power-law decay is preceded by a mean-field regime where the gap density decays with $`t^1`$. The two-stage relaxation and the mapping to the $`A+A0`$ model remain valid even if there is an attractive interaction between adsorbed dimers on the one-dimensional lattice. In the two-dimensional model with interaction, the sequential deposition stage goes over into a nucleation process and the relaxation stage into a much slower annealing process. Yet even then the two-stage relaxation remains qualitatively the same.

The annealing process for interacting $`k`$-mers with $`k>3`$ is fundamentally different from the dimer model. One finds that non-interacting $`k`$-mers are always well described with the mean-field kinetics of the model $`kA0`$ predicting a decay in the gap density $`t^{1/(k1)}`$. Introducing a strong interaction between the $`k`$-mers, as is the case for tropomyosin binding on actin filaments, we find a mapping to the $`A+AlA`$ model, valid for concentrations well above the equilibrium value, which shows a particle concentration decaying $`t^{1/2}`$. Closer to equilibrium other processes become important as well. First, particle splitting, $`AA+A`$ leads to an intermediate mean-field regime similar to that of non-interacting $`k`$-mers. Later the events of pair creation, $`0A+A`$ become relevant and lead the system into the final exponential relaxation towards equilibrium.

All of the above results concern the statics and dynamics of ligands decorating periodic structures in the absence of ATP hydrolysis. Such systems are “passive” in the sense that the ligands bind and unbind from their respective substrates but do not show any motor activity, i.e. move along these molecular tracks. As described above, this allowed us to give a detailed description of decoration experiments and to analyze corresponding experiments quantitatively. We also were able to identify and quantify cooperative effects between kinesin dimers. As our results show, the relaxation time in many experimental or biological systems can be extremely long. Therefore, one has to be aware that an experiment where the measurement is taken soon after the start of the decoration often does not show the equilibrium configuration, but rather some intermediate state from the relaxation process.

A natural extensions of this work is to investigate “active” systems composed of an ensemble of motor proteins and their respective substrates in the presence of transport along the molecular tracks driven by the chemical energy of ATP hydrolysis. Understanding such systems lies at the heart of many cellular processes. Although single processive motors like kinesin or myosin V can move loads over considerable distances, the number of motors acting simultaneously when transporting vesicles through the cytoplasm can be pretty large. For instance, in Ref. feneberg2001  the force acting on a bead in the cell has been estimated as $`100200\mathrm{pN}`$, which implies that at least $`2040`$ motor proteins were involved. How these active motors cooperate at high densities and what role is played by interactions between them, is to a large extend a still unexplored field with many open questions.

An even broader range of questions arises when one tries to understand how cells organize their interior to fulfill their various duties. While organizing its interior, the cell has to physically separate molecules or molecular aggregates from each other, define distinct functional regions, and actively transport molecules between these regions. Such processes rely on the assembly and ordering processes of macromolecules inside the cell and on forces that may be generated by molecular assembly or by the action of motor proteins. Understanding the principles underlying those complex cellular processes not only requires biological and biochemical studies but also studies of physical processes. In this respect microtubules and motors have been used as in vitro model systems to study various aspects of regulation and self-organization of cellular structures. It was shown theoretically and experimentally that the dynamic instability of microtubules in combination with microtubule polymerization forces is sufficient to provide a microtubule organizing center with a mechanism to position itself at the center of a confined geometry holy97 ; dogterom97 . In a simple in vitro system consisting solely of multi-headed constructs of the motor protein kinesin and stabilized microtubules the formation of asters and spiral defects was shown nedelec97 . By varying the relative concentration of the components a variety of self-organized structures (patterns) where observed. The pattern formation observed in these simple in vitro systems are only first examples of what will be a much broader range of self-tuned or regulated structural and temporal organization in biological systems. The control parameter determining the structure was the relative concentration of the components. In cellular processes there are, however, other possible parameters for controlling and regulating such as binding of a host of associated binding proteins. We expect that these and other related biological processes are a source for fascinating cooperative effects and a multitude of interesting nonequilibrium dynamic phenomena.

###### Acknowledgements.

We would like to thank Eckhard Mandelkow and Andreas Hoenger for helpful discussion about microtubule decoration experiments and their relevance for studying motor proteins. We have also benefited from discussions with Tom Duke, Thomas Franosch, Jaime Santos, Franz Schwabl, Gunter Schütz and Uwe C. Täuber. Our work has been supported by the DFG under contract nos. SFB 413 and FR 850/4-1. A.V. would like to acknowledge support by the European Union through a Marie Curie Fellowship (contract no. HPMFCT-2000-00522).


format: text

text:
# 1 Introduction

## 1 Introduction

Quantum field theory in curved space-time was actively developed in the 70-ties of the last century (see books?,?). However, there are a number of unsolved problems in this theory. One of them is related to the definition of an elementary particle and vacuum in a curved space-time. The reason is that in the case of a curved space there is no group of symmetries similar to the Poincaré group in the Minkowski space. If we believe that a particle is associated with a quantum of energy, then according to quantum mechanics, the observation of particles at an instant implies finding an eigenstate of the Hamiltonian. The Hamiltonian diagonalization method? takes this into account automatically.

The use of the Hamiltonian constructed from the metrical energy-momentum tensor, first proposed by A.A.Grib and S.G.Mamayev,? was successful in the homogeneous isotropic space-time for the conformal scalar fields. But such Hamiltonian leads to the difficulties related to an infinite density of created quasiparticles in the nonconformal case.? The energy of such quasiparticles differs from the oscillator frequency of the wave equation.?,? The nonconformal case is important for investigation by many reasons. In particular, the additional nonconformal terms are dominant in the vacuum expectation values of the energy-momentum tensor.? Recently,? the modified Hamiltonian was found that the density of the particles corresponding to its diagonal form and created in the nonstationary homogeneous isotropic space-time is finite.

The aim of this paper is the generalization of the method, proposed in Ref. 8, for the cases of the nongomogeneous and anisotropic spaces. In Sec. 2 we quantize the nonconformal scalar field in curved space and construct the modified Hamiltonian. In Sec. 3 we diagonalize this Hamiltonian and investigate the questions about the energies and the density of created particles.

## 2 Quantizing a Scalar Field in Curved Space

We consider a complex scalar field $`\phi (x)`$ of the mass $`m`$ with Lagrangian

$$L(x)=\sqrt{|g|}[g^{ik}_i\phi ^{}_k\phi (m^2+f(R))\phi ^{}\phi ],$$
(1)
and corresponding equation of motion

$$(^i_i+f(R)+m^2)\phi (x)=0,$$
(2)
where $`_i`$ are the covariant derivatives in $`N`$-dimensional space-time, $`g=\mathrm{det}(g_{ik})`$, $`f(R)`$ is the function of the invariant combinations $`g_{ik}`$ and curvature tensor $`R_{klm}^i`$. The equation (2) is conformally invariant if $`m=0`$ and $`f(R)=\xi _cR`$, where $`R`$ is scalar curvature, and $`\xi _c=(N2)/[\mathrm{\hspace{0.17em}4}(N1)]`$ (conformal coupling). We investigate the case of metric which in some coordinate system $`(x)=(\eta ,𝐱)`$ takes the form

$$ds^2=a^2(\eta )\left(d\eta ^2\gamma _{\alpha \beta }(x)dx^\alpha dx^\beta \right),$$
(3)
where $`\alpha ,\beta =1,\mathrm{},N1`$ and $`\mathrm{det}(\gamma _{\alpha \beta })=\gamma (𝐱)`$. It realizes, in particular, for the homogeneous isotropic space-time; for the nonhomogeneous spaces with

$$ds^2=dt^2a^2(t)\gamma _{\alpha \beta }(𝐱)dx^\alpha dx^\beta ;$$
(4)
for the anisotropic metric of Bianchi type-I $`\left(ds^2=dt^2a_\alpha ^2(t)(dx^\alpha )^2\right).`$ The equation (2) for the functions $`\stackrel{~}{\phi }=a^{(N2)/2}\phi `$ in coordinates $`(\eta ,𝐱)`$ takes the form

$$\stackrel{~}{\phi }^{\prime \prime }\mathrm{\Delta }_{N1}\stackrel{~}{\phi }+\left[\left(m^2+f(R)\right)a^2\frac{N2}{4}\left(2c^{}+(N2)c^2\right)\right]\stackrel{~}{\phi }=0,$$
(5)
where the prime denotes the derivative with respect to time $`\eta `$, $`ca^{}/a`$, $`\mathrm{\Delta }_{N1}\gamma ^{1/2}_\alpha (\sqrt{\gamma }\gamma ^{\alpha \beta }_\beta )`$. Let us assume that Laplace-Beltrami operator $`\mathrm{\Delta }_{N1}`$ has the complete orthonormal set of the eigenfunctions $`\mathrm{\Phi }_J(𝐱)`$ non-depending on time

$$\mathrm{\Delta }_{N1}\mathrm{\Phi }_J(𝐱)=\lambda ^2(J,\eta )\mathrm{\Phi }_J(𝐱).$$
(6)
The full set of the motion equation solutions can be found in the form $`\stackrel{~}{\phi }(x)=g_J(\eta )\mathrm{\Phi }_J(𝐱)`$ , where

$$g_J^{\prime \prime }(\eta )+\mathrm{\Omega }^2(\eta )g_J(\eta )=0,$$
(7)
$`\mathrm{\Omega }`$ is the oscillator frequency:

$$\mathrm{\Omega }^2(\eta )=\left(m^2+f(R)\right)a^2\frac{N2}{4}\left(2c^{}+(N2)c^2\right)+\lambda ^2(J,\eta ).$$
(8)
For quantization, we expand the field $`\stackrel{~}{\phi }(x)`$

$$\stackrel{~}{\phi }(x)=𝑑\mu (J)\left[\stackrel{~}{\phi }{}_{J}{}^{(+)}a_{J}^{(+)}+\stackrel{~}{\phi }{}_{J}{}^{()}a_{J}^{()}\right],$$
(9)
where $`d\mu (J)`$ is the measure on the space of $`\mathrm{\Delta }_{N1}`$ eigenvalues,

$$\stackrel{~}{\phi }{}_{J}{}^{(+)}(x)=\frac{1}{\sqrt{2}}g_J(\eta )\mathrm{\Phi }_J^{}(𝐱),\stackrel{~}{\phi }{}_{J}{}^{()}(x)=\left(\stackrel{~}{\phi }{}_{J}{}^{(+)}(𝐱)\right)^{},$$
(10)
and impose the usual commutation relations on the operators $`a_J^{(\pm )},\underset{J}{\overset{(\pm )}{\stackrel{}{a}}}`$.

We construct Hamiltonian as canonical one for variables $`\stackrel{~}{\phi }(x)`$ and $`\stackrel{~}{\phi }^{}(x)`$. If we add $`N`$-divergence $`(J^i/x^i)`$ to Lagrangian density (1), where in the coordinate system $`(\eta ,𝐱)`$ the $`N`$-vector $`(J^i)=(\sqrt{\gamma }c\stackrel{~}{\phi }^{}\stackrel{~}{\phi }(N2)/2,\mathrm{\hspace{0.17em}0},\mathrm{},0)`$, the motion equation (2) is invariant under this addition. By using the Lagrangian density $`L^\mathrm{\Delta }(x)=L(x)+(J^i/x^i)`$, we obtain by integration on hypersurface $`\mathrm{\Sigma }:\eta =const`$ of Hamiltonian density $`h(x)=\stackrel{~}{\phi }^{}(L^\mathrm{\Delta })/(\stackrel{~}{\phi }^{})+\stackrel{~}{\phi }^{}(L^\mathrm{\Delta })/(\stackrel{~}{\phi }^{})L^\mathrm{\Delta }(x)`$, the modified Hamiltonian

$`H(\eta )`$ $`=`$ $`{\displaystyle _\mathrm{\Sigma }}h(x)d^{N1}x={\displaystyle _\mathrm{\Sigma }}d^{N1}x\sqrt{\gamma }\{\stackrel{~}{\phi }^{}\stackrel{~}{\phi }{}_{}{}^{}+\gamma ^{\alpha \beta }_\alpha \stackrel{~}{\phi }^{}_\beta \stackrel{~}{\phi }+`$ (11)
$`+`$ $`[(m^2+f(R))a^2{\displaystyle \frac{N2}{4}}(2c^{}+(N2)c^2)]\stackrel{~}{\phi }^{}\stackrel{~}{\phi }\}.`$
Hamiltonian (11) is expressed in terms of the operators $`a_J^{(\pm )},\underset{J}{\overset{(\pm )}{\stackrel{}{a}}}`$ by

$$H(\eta )=d\mu (J)\{E_J(\eta )(\underset{J}{\overset{(+)}{\stackrel{}{a}}}a_J^{()}+\underset{\overline{J}}{\overset{()}{\stackrel{}{a}}}a_{\overline{J}}^{(+)})+F_J(\eta )\underset{J}{\overset{(+)}{\stackrel{}{a}}}a_{\overline{J}}^{(+)}+F_J^{}(\eta )\underset{\overline{J}}{\overset{()}{\stackrel{}{a}}}a_J^{()}\},$$
(12)
where

$$E_J(\eta )=\frac{1}{2}[|g_J^{}|^2+\mathrm{\Omega }^2|g_J|^2],F_J(\eta )=\frac{\vartheta _J}{2}\left[g_J^{}{}_{}{}^{2}+\mathrm{\Omega }^2g_J^{\mathrm{\hspace{0.17em}2}}\right],$$
(13)
here we used the normalization condition $`g_Jg_J^{}{}_{}{}^{}g_J^{}g_J^{}=2i`$ and the special choice of the eigenfunctions at which for an arbitrary $`J`$ such $`\overline{J}`$ existed, that $`\mathrm{\Phi }_J^{}(𝐱)=\vartheta _J\mathrm{\Phi }_{\overline{J}}(𝐱)`$.

## 3 Hamiltonian Diagonalization and Particle Creation

The Hamiltonian diagonalization for an arbitrary instant $`\eta `$ is realized in terms of the operators $`\underset{J}{\overset{(\pm )}{\stackrel{}{b}}},b_J^{(\pm )}`$ related to the operators $`\underset{J}{\overset{(\pm )}{\stackrel{}{a}}},a_J^{(\pm )}`$ via the time-dependent Bogolubov transformations

$$a_J^{()}=\alpha _J^{}(\eta )b_J^{()}(\eta )\beta _J(\eta )\vartheta _Jb_{\overline{J}}^{(+)}(\eta ),\underset{J}{\overset{()}{\stackrel{}{a}}}=\alpha _J^{}(\eta )\underset{J}{\overset{()}{\stackrel{}{b}}}(\eta )\beta _J(\eta )\vartheta _J\underset{\overline{J}}{\overset{(+)}{\stackrel{}{b}}}(\eta ),$$
(14)
where the functions $`\alpha _J,\beta _J`$ satisfy the initial conditions $`|\alpha _J(\eta _0)|=1,\beta _J(\eta _0)=0`$ and the identity $`|\alpha _J(\eta )|^2|\beta _J(\eta )|^2=1`$. Substituting expansion (14) in (12) and requiring that the coefficients at the nondiagonal terms $`\underset{J}{\overset{(\pm )}{\stackrel{}{b}}}b_J^{(\pm )}`$ vanish, we obtain

$$2\alpha _J\beta _J\vartheta _JE_J+\alpha _J^2F_J+\beta _J^2\vartheta _J^2F_J^{}=0,|\beta _J|^2=\frac{1}{4\mathrm{\Omega }}\left(|g_J^{}|^2+\mathrm{\Omega }^2|g_J|^2\right)\frac{1}{2},$$
(15)
$$H(\eta )=d\mu (J)\mathrm{\Omega }(\eta )(\underset{J}{\overset{(+)}{\stackrel{}{b}}}b_J^{()}+\underset{\overline{J}}{\overset{()}{\stackrel{}{b}}}b_{\overline{J}}^{(+)}).$$
(16)
Therefore, the energies of quasiparticles corresponding to the diagonal form of the Hamiltonian (12) are equal to the oscillator frequency $`\mathrm{\Omega }(\eta )`$.

The vacuum state for instant $`\eta `$ is defined by $`b_J^{()}|\mathrm{\hspace{0.17em}0}_\eta >=\underset{J}{\overset{()}{\stackrel{}{b}}}|\mathrm{\hspace{0.17em}0}_\eta >=0.`$ The state $`|\mathrm{\hspace{0.17em}0}>=|\mathrm{\hspace{0.17em}0}_{\eta _0}>`$ contains $`|\beta _J(\eta )|^2`$ quasiparticle pairs corresponding to the operators $`b_J^{(\pm )}(\eta )`$ in every mode.?,? The density of created particles is proportional to $`𝑑\mu (J)|\beta _J|^2`$. Using the first iteration in the integral equation for $`S(\eta )=|\beta _J(\eta )|^2`$

$$S(\eta )=\frac{1}{2}_{\eta _0}^\eta 𝑑\eta _1w(\eta _1)_{\eta _0}^{\eta _1}𝑑\eta _2w(\eta _2)(1+2S(\eta _2))\mathrm{cos}[2\mathrm{\Theta }(\eta _2,\eta _1)],$$
(17)
where $`w(\eta )=\mathrm{\Omega }^{}(\eta )/\mathrm{\Omega }(\eta )`$, $`\mathrm{\Theta }(\eta _1,\eta _2)=_{\eta _1}^{\eta _2}\mathrm{\Omega }(\eta )𝑑\eta `$, we can see, that if the eigenvalues $`\lambda ^2(J)`$ of the Laplace-Beltrami operator $`\mathrm{\Delta }_{N1}`$ don’t depend on time, then $`w\lambda ^2`$ and $`S\lambda ^6`$ in $`\lambda (J)\mathrm{}`$.

Therefore, in the space-time with metric (4) the density of created particles for $`N=4`$ is finite. (In Bianchi type-I the result is infinite). From equation (17) for $`S(\eta )`$ we can see that the densities of conformal and nonconformal (under $`f(R)=\xi R`$) particles created in Friedmann radiative-dominant Universe coincide.

Acknowledgements

The author is grateful to Prof. A.A.Grib for helpful discussions. This work was supported by Min. of Education of Russia, grant E00-3-163.

References