text:
# Comment on “Reexamination of experimental tests of the fluctuation theorem” by Narayan and Dhar

(July 11, 2003)

Our result in Ref.  that, instead of the conventional heat fluctuation theorem (FT), a new FT holds for heat fluctuations for a Brownian particle in a moving confining potential, was claimed to be disproved in a very recent preprint by Narayan and Dhar. This comment is meant to show that their assertion is not correct. The point is that they formulate their FT differently than we do ours. Effectively, their FT speaks about a physically irrelevant limiting case of our new FT. This implies that the two FTs are not in contradiction with each other. Furthermore, we point out an incorrect assumption in their derivation.

First, let us remark that although we have a FT for all times in Ref. , we only discuss the infinite time version in this comment.

What we consider as the conventional and the new FT, and what Narayan and Dhar in Ref.  consider the FT to be, respectively, can be formulated as follows.

We : $`P(Q_\tau =p\sigma \tau )`$ denotes the probability that the heat $`Q_\tau `$ produced in a time interval $`\tau `$, has the value $`p\sigma \tau `$, where $`\sigma `$ is the average heat produced per unit time ($`=Q_\tau /\tau `$) in the stationary state over positive times.***In Ref. , $`\sigma `$ is called $`w`$.

Narayan and Dhar : $`P(Q_\tau =A)`$ denotes the probability that the heat $`Q_\tau `$ produced in a time $`\tau `$ has the fixed value $`A`$, where $`A`$ is independent of $`\tau `$. We note that in Ref.  this is formulated in terms of $`S`$, which is the entropy production in a time $`\tau `$, and is given by $`S=Q_\tau /T`$, with $`T`$ the temperature of the fluid in which the Brownian particle moves.

As a consequence, there are two similarly looking FTs, but at closer inspection they are physically radically different.

We  consider the conventional FT in the form

$$\underset{\tau \mathrm{}}{lim}\frac{1}{\beta p\sigma \tau }\mathrm{log}\left[\frac{P(Q_\tau =p\sigma \tau )}{P(Q_\tau =p\sigma \tau )}\right]=1$$
(1)
for fixed $`p`$, where $`p`$ characterizes $`Q_\tau `$ in units of $`\sigma \tau =Q_\tau `$ and $`\beta =1/k_BT`$. This is consistent with Refs. .

Narayan and Dhar  consider the conventional FT to be

$$\underset{\tau \mathrm{}}{lim}\frac{1}{\beta A}\mathrm{log}\left[\frac{P(Q_\tau =A)}{P(Q_\tau =A)}\right]=1$$
(2)
for fixed, finite $`A`$.

Narayan and Dhar’s FT is so formulated — because of the independence of $`A`$ from $`\tau `$ — that when $`\tau `$ approaches infinity, so that $`Q_\tau `$ also approaches infinity, the fluctuation value $`A`$ remains always finite and therefore much smaller than $`Q_\tau `$. Then only values of $`A`$ much smaller than $`Q_\tau `$, but no values either near or larger than $`Q_\tau `$, are considered in Narayan and Dhar’s FT. Hence, only an increasingly restricted range of fluctuation values $`A`$ are actually studied in this FT as $`\tau `$ approaches infinity.

We, in contrast, consider all fluctuation values $`A`$, because we scale the value $`A`$ with $`\tau `$, so that when $`Q_\tau `$ grows ($`\tau `$), so does $`A=p\sigma \tau `$, as is written explicitly in our formulation in Eq. (1) for the FT. A $`\tau `$-proportional $`A`$ means that for a given $`Q_\tau `$, any heat fluctuations ($`Q_\tau =p\sigma \tau `$) can be considered by an appropriate choice of $`p`$, even as $`\tau `$ approaches infinity.This issue of scaling is analogous to using the proper scaling to derive, e.g., the diffusion equation for the random walk, which requires an $`l^2/t`$ scaling. Other scalings (like e.g., $`l/t`$) would not lead to a physically meaningful description of diffusion.

So we claim that Eq. (1) is consistent with Eq. (2) since a constant finite $`A`$ — as Narayan and Dhar’s FT uses in Eq. (2) — can also be obtained from Eq. (1), by taking the limit $`p0`$ (when $`\tau \mathrm{}`$). In other words, Eq. (2) is a very special case of Eq. (1) which covers strictly speaking only the value $`p=0`$, while Eq. (1) covers all values of $`p`$.

We will now argue that the FT of Narayan and Dhar in Eq. (2) is also consistent with our new FT found in Ref. . In Ref.  the conventional FT \[Eq. (1)\] was found to be incorrect as a FT for heat fluctuations in the system described above, and had to be replaced by a new FT. However, for the special region $`1p1`$, this new FT turns out to coincide with the conventional form of Eq. (1), which includes $`p=0`$. Therefore, there is no contradiction between our new FT and Narayan and Dhar’s FT in Eq. (2) either.

The above is enough to conclude that our new FT in Ref.  is not in contradiction with the result Eq. (2) (or their Eq. (1) in Ref. ). Moreover, our theory leading to the new FT in Ref.  has been compared in detail with two independent numerical calculations, one, a numerically carried out inverse Fourier transform explained in Ref. , and another, a sampling method, which is yet to be published. Both give perfect agreement with each other as well as with the new FT. So also from this, we have no reason to distrust the results in Ref. .

We take this opportunity to address very briefly two other points in Ref. .

Firstly, in their derivation of their (conventional) FT for $`Q_\tau `$ in their Sec. III, Narayan and Dhar neglect a term in their Eq. (17) which one can neglect for the average $`Q_\tau `$, but not for the fluctuations of $`Q_\tau `$. In fact, it is precisely this term which gives rise to exponential rather than Gaussian behavior of the tails of the distribution $`P(Q_\tau )`$ — the existence of which Narayan and Dhar also notice — and through these exponential tails, this term also gives rise to our new FT, instead of to the conventional one. To substantiate this, a saddle point method, involving a complicated study of the behavior of the Fourier transform of $`P(Q_\tau )`$, in particular, its singularities in the complex plane, is needed, which goes beyond the scope of this comment.

Secondly, Narayan and Dhar seem to suggest that the singularity of the Langevin equation for vanishing mass $`m`$ of the particle might be a problem in the derivation of a FT. However, setting $`m=0`$ from the start in their Eq. (12) or setting $`m=0`$ at the end in their Eqs. (17) and (19), give the same FT, i.e., no singularity is encountered in this derivation with respect to $`m=0`$.


format: text

text:
# Simulations of Oligomeric Intermediates in Prion Diseases

## I INTRODUCTION

Prion diseases are a group of neurodegenerative diseases including bovine spongiform encephalopathy (BSE) in cattle, scrapie in sheep and goats, chronic wasting disease in deer and elk, and kuru and Creutzfeldt-Jakob Disease (CJD) in humans. These diseases came to the forefront after BSE reached epidemic proportions in Great Britain in the early 1990s, and it was later shown that transmission of BSE to humans can lead to new variant CJD (vCJD) in humans (Bruce et al., 1997), (Hill et al., 1997), (Scott et al., 1999).

Prion diseases are unusual in that they appear to be caused by infection with some minimal infectious ”seed” of misfolded prion protein, which alone may be able to cause disease by catalyzing further misfolding and, in many cases, aggregation of the prion protein. These aggregates are typically amyloid-like fibrils or amyloid plaques (Caughey, 2000). The infectious agent is unusually hard to eliminate by various methods including UV irradiation, suggesting it contains no nucleic acid and ratheronly protein, the so-called ”protein-only” hypothesis in prion diseases(Weissmann et al., 2002).

In the case of CJD, a sporadic form of the diseases also exists, occurring more or less randomly worldwide with an incidence of about one in a million people per year. It has been suggested that this incidence is due to the very rare event of nucleating the minimal infectious seed by chance in a healthy individual (Come et al., 1993).

Developing an understanding of these diseases is important because, for one, they are invariably fatal. To date, no treatment exists. Additionally, it is not yet clear how large the vCJD epidemic in humans will be; an understanding of the disease process is important to be able to guide the search for treatment ideas.

In many cases, prion diseases result in large, up to micron-scale plaques in the brains of people and animals with these diseases. They also involve vacuolization or spongiform change in the brain due to death of neurons(Scott et al., 1996). Additionally, the normal form of the prion protein (known as PrPC) has long been known to misfold and aggregate *in vitro* when catalyzed by the presence of a misfolded prion protein (PrPSc) seed (Come et al., 1993). Together, these observations have suggested to some that the aggregation process itself may be important in these diseases (Come et al., 1993; Masel et al., 1999). It has also been suggested that the rate-limiting step in aggregation is nucleation of an appropriate seed, thus the rapid aggregation in the seeded case described above (Come et al., 1993).

Another fact which may be important to this issue is that the prion protein is normally GPI-anchored to the cell surface. Aggregation  *in vitro* as mentioned above is observed in solution rather than in the presence of the GPI anchor on a cell surface, leaving the possibility that the aggregation process in vivo is different.

Aggregation models developed to explore the aggregation process in prion disease include one-dimensional, fibrillar aggregation-and-fission models(Masel et al., 1999; Slepoy et al., 2001), since aggregates grown *in vitro* are typically seen to be fibrillar. Additionally, our earlier work suggested that a two-dimensional areal aggregation model could explain certain other properties of the diseases (Slepoy et al., 2001). This earlier model is attractive in that it can provide a simple explanation for the long lag phase which is sometimes observed in growth of the amount of infectious material in the brain. This lag phase of little or no growth is followed by a doubling phase with a short characteristic doubling time. Additionally, our earlier model provides a possible explanation of some of the difference between infectious and sporadic forms of CJD (Slepoy et al., 2001). In later work, we used this model to explain and fit experimental dose incubation curves (Kulkarni et al., 2002).

However, there were drawbacks to the earlier aggregation model we proposed. First, no such areal aggregates had so far been observed. Second, the fissioning essential to the model would involve breaking of strong bonds between the proteins, probably bonds between $`\beta `$-sheets (Serag et al., 2002).

More recent experimental work found two-dimensional areal aggregates of prion protein produced during the purification process. These aggregates were examined under electron microscope (EM) and found to consist trimeric or hexameric subunits. These subunits are linked together in a regular array, possibly by their N-terminal sugars or a weak protein-protein interaction (Wille et al., 2002).

This suggested we should modify our earlier model and attempt to reproduce this aggregate morphology. We thought of two basic schemes for growing aggregates of this sort: (1) Growing the aggregate outward monomer by monomer from an initial seed, or (2) Oligomeric intermediates (possibly very flexible and of unstable shape), form on their own in solution and are only catalyzed into stably misfolding in the presence of an existing misfolded seed. Some evidence in favor of (2) has already been produced: Monomers of yeast prion can form intermediates if left to stand, which allows aggregation to proceed at an initial faster rate when catalyzed by addition of a seed (Serio et al., 2000). Additionally, the conformation-dependent immunoassay developed by Safar et al. (2002) detects both protease-sensitive and protease-resistant PrPSc. In hamster brains, sensitive PrPSc is observed earlier, followed by resistant PrPSc. This could correspond to (2) above, where the sensitive PrPSc is the intermediates that are not yet stably misfolded and the resistant PrPSc is stably misfolded intermediates.

Work here has been done to further explore these two potential modifications of our earlier model to examine whether they retain the same features and if additional insight can be gained.

It is important to note that even if areal aggregation is not important to the time course of these diseases, the aggregates observed by Wille et al. have already provided insight into the structure of the misfolded prion protein (Wille et al., 2002). Theoretical modeling may be able to place further constraints on the protein or subunit structure necessary to reproduce these aggregates, and hence provide valuable information because these aggregates *can* form, even if they are not important to the disease progression.

## II BASICS OF OUR MODEL

Here we explore the two basic schemes suggested above for growing aggregates like those observed by Wille et. al. To do so, we use a modification of our earlier model. Therefore a recap of common features of these models is useful.

These models are stochastic cellular automata models, meaning that they take place on a lattice with probabilistic interaction and diffusion rules governing the progression of the system. In this case, sites on the lattice are either occupied by individual prion proteins, or water (empty, in the simulation). The protein form at a site can also vary from PrPC to PrPSc.

Rules vary depending on the model being explored, but the basic procedure is the same. Every simulation step, which represents a small amount of time, we allow proteins and any aggregates to diffuse a small amount on the lattice (each object has a probability $`1/(size)^{1/2}`$ of moving one lattice site in a given step). Then we look at every protein in the lattice and update its state according to the rules. For example, in our original model, the conformation of an individual prion protein is determined solely by its number of neighboring prion proteins, and this can vary from step to step. After doing this, we add more normal prion monomers to replace any that converted to PrPSc. This is due to the assumption that this process would be taking place in a small area on a cell, and the normal prion monomers would be added by the cell or diffuse in from other locations on the cell surface to keep the monomer concentration relatively constant.

## III MONOMER GROWTH MODEL

First, case (1) from above was explored. Simple rules were developed (figure 1) which can reproduce aggregates similar to those observed by Wille et al. It is important to note that though the rules were designed to reproduce such aggregates, most modifications of these rules could not do so. This means that the rules provide some constraints on the protein-protein interactions necessary to reproduce such aggregates. Also, for the purposes of this model, we are assuming the subunits are hexameric, but the corresponding model for trimeric intermediates is actually much simpler than this model and will produce similar results.

Details of the algorithm for this model are covered in figure 2.

The rules are as follows: The simulation begins with a single hexagonal subunit consisting of six misfolded monomers (light gray hexagons in figure 1) which stick some of their residues into an adjacent site, excluding anything else from occupying that site (black). Healthy monomers (light gray spheres) can then attach via a sugar-bond or other weak protein-protein interaction to this subunit (dark gray spheres to dark gray hexagons) but only radially outward from a monomer in the initial hexamer. Additional monomers moving adjacent to the attached monomer can, together with it, misfold but only if the second monomer does not also neighbor the original hexamer. Then additional monomers can attach to this forming hexamer, allowing it to complete. Repeating this process many times can produce mostly regular aggregates with some holes, similar to those observed.

The rules are also probabilistic: above, ”can” means that some fraction of the time the event occurs. These probabilities can be changed in the simulation and give different growth rates, but the same essential features and scaling as described below.

If this is in fact how these aggregates are forming, we find out about the orientation of monomers within a hexagonal subunit. We find, as mentioned in the discussion of the rules above, that the N-terminal sugars or attachment sites must stick radially outward from each monomer in a hexagonalsubunit (figure 1b). This is in agreement with the hexagonal structure proposed by Wille. Additionally, we find that no such regular aggregates can be produced unless the monomer attaching to a previously attached monomer (figure 1(c)) can only attach if it is not adjacent to an existing hexamer. This seems to indicate that the other spaces must be occupied by residues from the existing hexamer, preventing attachment in those sites.

This model can also reproduce gaps in aggregates as observed. In this model gaps are due to variations of the growth rate from average for part of the aggregate, causing several parts of the aggregate to grow apart and then rejoin after leaving a gap.

One reason for developing this model was to see if it would capture the same features of the disease as our original model. Our original model explained the difference between the lag phase and the doubling phase by suggesting that the doubling phase is initiated when aggregates begin to fission, then regrow to a certain fissioning size and break again. Key to this explanation is our result that aggregation speeds up, so that the time for an aggregate to double in size from half its fission size to its fission size is much much less than the time for it to get from its initial size to its fissioning size.

To see if this model could produce the same separation of lag and doubling phases, we examined the aggregate growth rate as a function of size in this model (figure 3) and found it speeds up only slowly. Naïvely, one would expect the growth rate to be roughly proportional to the square root of the size, as the growth rate is proportional to the circumference of the aggregate, which, assuming a circular aggregate, is $`2\pi r`$. The size of the aggregate is proportional to the area, $`\pi r^2`$, so the radius is proportional to the square root of the size and thus the rate proportional to the square root of the size. To a good approximation, the growth rate observed here is well-fit by an offset plus a term proportional to $`(size)^{1/2}`$, as expected.

In this simple picture, one can calculate the ratio of the doubling time to the lag time. The lag time is the time to go from the initial size, say size 0 for simplicity, to size n; the doubling time from size n/2 to size n. Integrating the rate to get the times and taking the ratio we find $`t_{doub}/t_{lag}=11/\sqrt{2}`$ or about 0.293. This means that this model cannot produce such a large separation between lag and doubling times as our earlier model could, at least without further modification.

This also indicates that if there is a lag phase and if the difference between it and the doubling phase is due to acceleration of aggregation, this picture is not sufficient and something more like (2), growth from intermediates, may be a better representation of the disease process.

## IV GROWTH VIA INTERMEDIATES

In this case, aggregation is assumed to be the assembly of independent hexameric intermediates into a larger areal aggregate. The intermediates themselves are not misfolded but only misfold, in this model, when they either aggregate with an existing misfolded seed, or come together in such a way that they can misfold and and form a new stable seed. In this way, the model works essentially just like the model of Slepoy et al. except now hexameric intermediates are playing the role of monomers (figure 4). As mentioned above, there is some evidence that intermediates greatly increase aggregation rate in studies of yeast prions, so this emphasis on the importance of intermediates may be reasonable.

To be able to map this model back into our old model, though, it must be known how the intermediate concentration depends on monomer concentration. And this is not obvious. So a simulation was developed to explore how the concentration of hypothetical hexameric intermediates would depend on monomer concentration. Again, here we are assuming the intermediates are hexameric but we can easily modify the model to accomodate trimers.

To get at the concentration of intermediates, it was assumed that two monomers have a probability $`P_1`$ of beginning a new hexameric subunit when they come into contact (see figure 5). This new subunit can grow by addition of monomers when they move into appropriate positions (changing this probability does not affect the outcome of the simulation, only the timescale, so it was set to 1). However, this growth process competes with a ”dissolving” process by which a monomer that is part of an intermediate but only has one neighboring monomer can break off with a probability $`P_3`$. Thus the end destiny of any intermediate that begins is either to form a complete hexameric intermediate, in which case it can persist, or to dissolve completely.

Details of the algorithm for this model are shown in figure 6.

This dissolving, or reversibility, was included because it was not obvious that at low monomer concentrations, one would expect a reasonable formation rate of intermediates via this mechanism. It was initially thought that at concentrations below something on the order of $`P_3`$, breaking would dominate and the formation rate of intermediates would be almost zero. First, the simulation that was developed was used to examine the dependence of time for intermediate formation as a function of monomer concentration (figure 7). It was found that at high monomer concentration, the time to form an intermediate scales between $`1/c`$ and $`1/c^2`$ ($`c`$ is concentration). This is because the likelihood of starting an intermediate scales as the dimer concentration ($`1/c^2`$) while the time to add monomers to it scales as $`1/c`$. On the other hand, at very low monomer concentration, the time asymptotically approaches $`1/c^6`$. This is due to the fact that at these concentrations, dissolving dominates and it is only in the very rare event that six monomers are in the same place at almost the same time that an intermediate can finish. The probability of that scales as $`1/c^6`$.

It is interesting to note that the beginning of the transition between high concentration behavior, where most intermediates successfully become complete, and low concentration behavior, where only a lucky few do, begins at a concentration on the order of the breaking probability, $`P_3`$. This suggests that if the strength of bonds between intermediates could be weakened somehow, the biological number of intermediates could be drastically decreased by pushing biological monomer concentrations into the $`1/c^6`$ regime.

The goal, however, was to determine the dependence of the intermediate concentration on monomer concentration. This just provided a formation rate, and the functional form was uncertain. So another sort of result was examined: We began examining behavior of the system as a function of time, and measured the number of different partial intermediates (two monomers, …five monomers, hexameric intermediates). We first examined the case with no breaking ($`P_3=0`$) to check our results, because it is relatively easy to work out kinetics in that case. A sample of one of these plots is shown in figure 8, with symbols as data points and solid lines as approximate kinetics fits. It is important to note that in this case, and in the case of nonzero breaking probability, the number of dimers, trimers, tetramers, and pentamers reaches equilibrium relatively quickly and then the hexamer number begins to grow linearly at a rate equal to the rate of dimer formation.

Sample results with nonzero breaking are shown in figure 9. These results are qualitatively similar, except the number of pre-intermediates that persists is much lower. In the high-breaking limit, the very low level of intermediates demonstrates that either a potential intermediate gets “lucky” and quickly forms an intermediate, or it dissolves back to monomers, leaving few dimers, trimers, and so on.

The kinetics equations we can write down to describe this simulation are relatively simple. With $`r_{nm}`$ the rate constant for forming m-mers from n-mers, and $`b_{nm}`$ the rate of breaking n-mers into m-mers plus monomers, we can write:

$$[c_1]=c$$
(1)
$$\frac{d[c_2]}{dt}=r_{12}[c_1]^2r_{23}[c_1][c_2]b_{21}[c_2]+b_{32}[c_3]$$
(2)
$$\frac{d[c_3]}{dt}=r_{23}[c_1][c_2]r_{34}[c_1][c_3]b_{32}[c_3]+b_{43}[c_4]$$
(3)
$$\frac{d[c_4]}{dt}=r_{34}[c_1][c_3]r_{45}[c_1][c_4]b_{43}[c_4]+b_{54}[c_5]$$
(4)
$$\frac{d[c_5]}{dt}=r_{45}[c_1][c_4]r_{56}[c_1][c_5]b_{54}[c_5]$$
(5)
$$\frac{d[c_6]}{dt}=r_{56}[c_1][c_5]$$
(6)
Since we know that the hexamer number grows linearly at steady state and all of the other concentrations are unchanging, we can greatly simplify the above kinetics by looking at the steady state only. We can work backwards from the steady state behavior of the hexamers to find the dependence of the steady state rate of hexamer formation on the different kinetic parameters and ultimately on the monomer concentration.

This straightforward kinetics analysis produces the equlibrium result

$$m=\frac{r_{12}c^2}{1+\frac{b_{21}}{r_{23}c}\left\{1+\frac{b_{32}}{r_{34}c}\left[1+\frac{b_{43}}{r_{45}c}\left(1+\frac{b_{54}}{r_{56}c}\right)\right]\right\}}$$
(7)
where $`m`$ is the slope at equilibrium of the hexamer formation rate.

The constants in our simple result for $`m`$, above, can be measured from our simulation. However, our simulation does not necessarily reproduce what these constants would be in a biological system. So it is difficult to say exactly what the rate of intermediate formation, $`m`$, would be in a real system. However, it is nevertheless useful to know the functional form of its dependence on the monomer concentration.

The result that the hexamer number begins growing linearly eventually is independent of monomer concentration. This is important because it means some hexamers can form given these simple rules even if breaking dominates. Given that result, it seems safe to assume that if hexameric intermediates are stable, some will form in biological systems.

In our model, the hexamer number grows linearly indefinitely, which is obviously unrealistic biologically. The reason for this is that we include no mechanism to remove finished hexamers. Realistically, they would be cleared from the body somehow. They could be endocytosed from the cell surface and degraded via the proteasome mechanism or some other pathway. Additionally, any hexamers being taken up into aggregates would reduce this number. Regardless, realistically the number should stabilize at some fixed value determined by the balance of the clearance rate and the formation rate.

With the result that some hexamers form even at low monomer concentrations (and more would form if they are trimers), a model was developed where now hexameric intermediates occupy a single cell on the lattice (equivalently, these could be trimeric intermediates). This model, described below, largely maintains the same attractive features of the original, showing that if areal aggregation is the explanation for these features, as we suggested, this aggregation could be of hexameric intermediates.

Part of our basis for this model is the observation that the intermediates are not yet stably misfolded since formation of intermediates in studies of yeast prions does not lead to a change in circular dichroism results; it is only when they aggregate with a seed that they stably misfold (Serio et al., 2002). If any aggregates consist of oligomers like those observed by Wille et al., then any misfolding of these intermediates must be catalyzed by existing aggregates or few-hexamer oligomers. We hypothesize that the mechanism for this is intermediates forming bonds to an existing seed. When solvent is excluded locally around these oligomers and their neighbors include a misfolded oligomer or aggregate, they misfold. The important point is that it is solvent exclusion around an intermediate that can cause it to misfold, making this a very rare sporadic event. But a misfolded seed can help this process by providing a place where intermediates bond, helping the solvent exclusion process. These rules make this model essentially identical in terms of kinetics to our original model.

Details of the algorithm for this model and mapping are shown in figure 10.

However, from our old model we estimated the sporadic form of the disease could have a peak around 1000 years given a biological concentration of $`10^3`$%. In our new model we find that it is very difficult to estimate this number as the scaling of the time as a function of monomer concentration is complicated. It was hoped that this model would give a result for the onset of sporadic disease that could be compared with the time for onset of the infectious form to see if the results were consistent with the roughly 1 in $`10^6`$ incidence of sporadic CJD that we earlier pointed out. Unfortunately, it is difficult for our model to give a concrete answer at this time as the answer depends too much on the value of the biological monomer concentration. We do find, however, that the power law used previously to scale the sporadic data, $`c^3`$, is a lower bound on the separation. That is, the actual exponent should be larger, meaning that we previously underestimated the separation of timescales. Thus although we cannot say exactly what the separation of timescales here will be, we can say that it will be greater than the two orders of magnitude that we previously estimated.

This work suggests that a model like our earlier one, modified to involve areal aggregation of hexameric or trimeric intermediates, could maintain the same attractive features of our earlier model in explaining certain aspects of the diseases. However, without precise knowledge of the biological monomer concentration and a way to measure relevant rate constants, it is difficult to make numerical predictions from this model.

## V DISCUSSION

Our work has shown that both in the case of monomer addition to a seed, and in the case of growth via intermediates, it is possible to produce aggregates like those observed by Wille et al. This leaves the question of how such aggregates actually grew. If areal aggregation is the cause, or part of the cause, of the difference between lag and doubling times, as suggested by Slepoy et al., then our work suggests that intermediates are already present *in vivo* prior to aggregation.

Our work has also shown that a model can be developed which, with suitable parameters, can reproduce areal aggregates like those actually observed while maintaining the same features of our original model.

Whether or not areal aggregation is actually important in these diseases, we can gain insight from this model. If the aggregates observed are growing via monomer addition, we gain some constraints on the structure simply from our rules. On the other hand, if intermediates are important to aggregation, then our results indicate the intermediate concentration can be quite important. At high intermediate concentrations, intermediates form relatively fast. However, at low intermediate concentrations, intermediate formation time scales as $`1/c^6`$. This result is exciting because it suggests the intermediates as a target to prevent aggregation. Simply reducing the monomer concentration by a factor of 2 would decrease the number of intermediates by a factor of $`2^6`$ or 64. Within our model, this would certainly increase the aggregation time, and thus slow down the disease, by at least the same factor. For a disease which can typically incubate for years, this obviously would be a great advantage.

In this case, the location of the transition between low concentration behavior and high concentration behavior is, roughly speaking, set by the probability of monomers breaking off from an intermediate before it becomes a stable hexamer. Thus if this probability could be increased slightly — that is, the bonds between monomers could be weakened slightly — it would have the result described above. This could provide an explanation for one experimental observation. Humans have a methionine/valine polymorphism at codon 129 of the gene for the prion protein. To date, everyone affected by vCJD has been methionine/methionine homozygous. This effect was also seen in the prion disease Kuru, where the methionine/methionine genotype was associated with increased susceptibility and the shortest incubation time (Goldfarb, 2002). If replacing methionine with valine weakened the monomer-monomer bonds within a forming intermediate and reduced intermediate concentration, this could have exactly the effect described above.

In all, our work shows that our earlier model can be extended to produce aggregates like those observed *in vitro* while still maintaining its attractive features. Our work also suggests possible mechanisms for formation of these aggregates. If the aggregates form by monomer addition, it constrains protein structure. If they form by addition of intermediates, it highlights the importance of bonds within the intermediates as a target for possible treatment strategies. Our model suggests that an experiment to measure the biological intermediate concentration, if there is such a concentration, would be very useful. That would indicate whether such intermediates are present at a high enough concentration to be important biologically. Additionally, this work suggests that experimentalists should check and see whether reasonably-sized aggregates of prion protein can be found *in vivo* on the cell surface. This confinement to the cell surface concievably could make the difference between the 1D fibrillar aggregates typically observed *in vitro* and 2D areal aggregates like those suggested by the model of Slepoy et al. Direct measurements, or detailed simulations, giving the strengths of beta bonds between monomers compared to bonds between subunits would be very useful.

One general scheme for experimentally testing the possible role of intermediates and estimating their concentration is via spin labelling (Hubbell 1998, Columbus 2002). Briefly, a small molecule with a free spin can preferentially react and attach to cysteine residues in proteins. Frequently, these residues are moved around a protein via mutagenesis to then map out structures, but for these purposes a less refined approach is required. Since the PrP protein already possesses cysteine residues at the position of the disulfide bond, the spin labels can attach there (and will not disrupt the disulfide bond). Then the spin-spin interactions will produce a different characteristic spectrum for monomers, incomplete intermediates and complete intermediates, in particular, with a progressive broadening upon moving from monomers to complete intermediates. Since the spins can have interactions with other spins within a 3 nm sphere, we don’t doubt that the broadening will be observable. Of course, since the spin labels will react with *any* cysteines present, it is important to carry this out first in *in vitro* aggregation experiments with purified prion extracts. This will help to identify the conditions which can lead to areal aggregation as observed by Wille et al. (2000), and serve as an existence proof at least for significant oligomeric intermediate concentrations.

This research is supported in part by the NEAT-IGERT program sponsored by the National Science Foundation (IGERT Grant DGE-9972741).

We gratefully acknowledge fruitful discussions on spin labelling with John Voss, and on areal prion aggregates with Holger Wille. R.R.P.S. and D.L.C. have benefitted from discussions at workshops of the Institute for Complex Adaptive Matter.

Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy’s National Nuclear Security Administration under contract DE-AC04-94AL85000.

References:

Bruce, M.E., Will, R.G., Ironside, J.W., McConnell, I., Drummond, D., Suttie, A., McCardle, L., Chree, A., Hope, J. Birkett, C., Cousens, S., Frasier, H., Bostock, C.J. 1997. Transmissions to mice indicate that ’new variant’ CJD is caused by the BSE agent. *Nature* 389: 498-501.

Caughey, B. 2000. Prion protein interconversions. *Philos. Trans. R. Soc. Lond. B Biol. Sci.* 356:197-202.

Columbus, L. and Hubbell, W.L. 2002. A new spin on protein dynamics. *Trends Biochem. Sci.* 27: 288-295.

Come, J.H., Fraser, P.E., Lansbury, P.T. Jr. 1993. A kinetic model for amyloid formation in the prion diseases: Importance of seeding. *Proc. Nat. Acad. Sci. USA* 90: 5959-5963.

Goldfarb, L.G. 2002. Kuru: The old epidemic in a new mirror. *Microbes Infect.* 4: 875-882.

Hubbell, W.L., Gross, A., Langen, R., and Lietzow, M.A. 1998. Recent advances in site-directed spin labelling of proteins. *Curr. Opin. Struct. Biol.*, 8: 649-656.

Hill, A.F., Desbruslais, M., Joiner, S., Sidle, K.C.L., Gowland, I., Collinge, J., Doey, L.J., Lantos, P. 1997. The same prion strain causes vCJD and BSE. *Nature* 389: 448-450.

Kulkarni, R., Slepoy, A., Singh, R.R.P., Pazmandi, F. 2003. Theoretical modeling of prion disease incubation. *Biophys. J.* In press.

Masel, J., Jansen, V.A.A., Nowak, M.A. 1999. Quantifying the kinetic parameters of prion replication. *Biophys. Chem.* 77: 139-152.

Philips, Lord of Worth Malraven, Bridgeman, J., Ferguson-Smith, M. October 2000. *Report of the BSE Inquiry*, Volume 1, Chapter 12, Paragraph 1122(i). http://www.bseinquiry.gov. Jan. 4, 2003.

Safar, J.G., Scott, M., Monaghan, J., Deering, C., Didorenko, S., Vergara, J., Ball, H., Legname, G., Leclerc, E., Solforosi, L., Serban, H., Groth, D., Burton, D.R., Prusiner, S.B., Williamson, R.A. 2002. Measuring prions causing bovine spongiform encephalopathy or chronic wasting disease by immunoassays and transgenic mice. *Nature Biotech.* 20: 1147-1150.

Scott, M.R.D., Telling, G.C., Prusiner, S.B. 1996. Transgenetics and Gene Targeting in Studies of Prion Diseases. *In* Prions Prions Prions, S.B. Prusiner, editor. Springer-Verlag, Berlin, Heidelberg. 95-123.

Scott, M.R., Will, R., Nguyen, H.-O. B., Tremblay, P., DeArmond, S., Prusiner, S.B. 1999. Compelling transgenetic evidence for transmission of bovine spongiform encephalopathy prions to humans. *Proc. Nat. Acad. Sci. USA* 96(26):15137-15142

Slepoy, A., Singh, R.R.P, Pázmándi, F., Kulkarni, R.V., Cox, D.L. 2001. Statistical Mechanics of Prion Diseases. *Phys. Rev. Lett.* 87(5): 581011 to 581014.

Serag, A.A., Altenbach, C., Gingery, M., Hubbell, W.L., Yeates, T.O. 2002. Arrangement of subunits and ordering of beta-strands in an amyloid sheet. *Nature Struct. Biol.* 9: 734-739.

Serio, T.R., Cashikar, A.G., Kowal, A.S., Sawicki, G.J., Moslehi, J.J. Serpell, L., Arnsdorf, M.F., Lindquist, S.L. 2000. Nucleated conformational conversion and the replication of conformational information by a prion determinant. *Science* 289: 1317-1321.

Wille, H., Michelitsch, M.D., Guenebaut, V., Supattapone, S., Serban, A., Cohen, F.E., Agard, D.A., Prusiner, S.B. 2002. Structural studies of the scrapie prion protein by electron crystallography. *Proc. Nat. Acad. Sci. USA* 99(6): 3563-3568.

Weissmann, C., Enari, M., Klohn, P-C., Rossi, D., Flechsig, E. 2002. Transmission of Prions. *J. Infect. Dis.* 186(Suppl. 2): S157-65.


text:
# Chiral Dynamics in the Meson Sector at two Loops

## 1 Introduction

The problem of dynamics of mesons at low energies is important. It plays a major role in the precise determination of the elements of the Cabibbo-Kobayashi-Maskawa matrix (CKM) which is a main part of the study of the standard model flavour sector CKM . In this talk I will concentrate on the theory behind the measurement of $`V_{us}`$ from $`K_\mathrm{}3`$ ($`K\pi \mathrm{}\nu `$) decays and in particular on the recent work of P. Talavera and myself on the two-loop calculation and the determination of the relevant low-energy constants BT2 ; Moriond .

I review shortly Chiral Perturbation Theory (ChPT) and the relevant two-loop calculations done up to now, followed by a discussion of $`K_\mathrm{}3`$ and the present theory situation. I include a short discussion of the validity of the linear approximation of the form factors normally used in the data analysis. The main relevant results from ChPT BT2  can be summarized as follows. The curvatures are important in the analysis but can be predicted using ChPT from the pion electromagnetic form-factor BT2  and all order $`p^6`$ parameters needed to determine $`V_{us}`$ can be experimentally obtained via the scalar form factor, $`f_0(t)`$, in $`K_\mathrm{}3`$ BT2 . The curvature of $`f_0(t)`$ can be predicted as well from knowledge about scalar form factors of the pion BD , albeit only at fairly low precision at present.

## 2 Chiral Perturbation Theory

ChPT is an effective field theory valid as an approximation to Quantum Chromodynamics (QCD) at low energies. Its modern form has was introduced by Weinberg, Gasser and Leutwyler ChPT ; GL1 . The global chiral symmetry, $`SU(3)_L\times SU(3)_R`$, of QCD in the limit of massless quarks is spontaneously broken down to the diagonal subgroup $`SU(3)_V`$ by a nonzero quark condensate, $`\overline{q}q=\overline{q}_Lq_R+\overline{q}_Rq_L0.`$ The eight broken generators lead to eight Goldstone bosons. These are massless and their interactions vanish at zero momentum. The latter allows the construction of a well defined perturbative expansion in terms of momenta, generically referred to as an expansion in $`p^2`$. Quark masses are usually counted as order $`p^2`$ since $`p_\pi ^2=m_\pi ^2m_q\overline{q}q/F_\pi ^2`$. Inserting an external photon or $`W^\pm `$-bosons counts as order $`p`$ since these are included via covariant derivatives. Recent lectures, much more detailed than what is included here are chptlectures . ChPT being an effective field theory implies that the number of parameters increases order by order. In the purely mesonic strong and semi-leptonic sector there are two parameters at lowest order ($`p^2`$), ten at NLO ($`p^4`$GL1 , and 90 at NNLO ($`p^6`$) BCE1 . The renormalization procedure and the divergences are worked out in general to NNLO BCE2  and provide a good check on all calculations. One problem shared with other high order loop calculations in comparing different calculations is the use of different renormalization schemes. The calculations that were used to determine all the needed parameters are those of the masses and decay constants ABT1 , $`K_\mathrm{}4`$ ABT2  and the electromagnetic form factors BT1 .

## 3 $`K_\mathrm{}3`$: definitions, $`V_{us}`$ and form factor linearity

The neutral and charged $`K_\mathrm{}3`$ decays amplitudes, $`K^{+,0}(p)\pi ^{0,}(p^{})\mathrm{}^+(p_{\mathrm{}})\nu _{\mathrm{}}(p_\nu )`$, are

$`T^{(+,0)}`$ $`=`$ $`{\displaystyle \frac{G_F}{\sqrt{2}}}V_{us}^{}\mathrm{}^\mu F_{\mu }^{}{}_{}{}^{(+,0)}(p^{},p),\mathrm{}^\mu =\overline{u}(p_\nu )\gamma ^\mu (1\gamma _5)v(p_{\mathrm{}}),`$
$`F_{\mu }^{}{}_{}{}^{+,0}(p^{},p)`$ $`=`$ $`(1/\sqrt{2},1)[(p^{}+p)_\mu f_+^{K^+\pi ^0,K^0\pi ^{}}(t)+(pp^{})_\mu f_{}^{K^+\pi ^0,K^0\pi ^{}}(t)].`$ (1)
Isospin leads to the relations

$$f_+^{K^0\pi ^{}}(t)=f_+^{K^+\pi ^0}(t)=f_+(t)\text{and}f_{}^{K^0\pi ^{}}(t)=f_{}^{K^+\pi ^0}(t)=f_{}(t),$$
(2)
The scalar form factor and the usual linear parameterizations are defined as

$$f_0(t)=f_+(t)+t/(m_K^2m_\pi ^2)f_{}(t),f_{+,0}(t)=f_+(0)\left(1+\lambda _{+,0}t/m_\pi ^2\right).$$
(3)
To determine $`|V_{us}|`$ we need $`f_+(0)`$ theoretically and experimentally. There are three main theoretical effects. There is a well-known short-distance correction from $`G_\mu `$ to $`G_F`$ calculated by Marciano and Sirlin. The corrections of order $`(m_s\widehat{m})^2`$ allowed by the Behrends-Sirlin-Ademollo-Gatto theorem are discussed at $`p^6`$ here and in BT2 ; Moriond . The sizable isospin breaking found by Leutwyler and Roos LR  is in the process of being evaluated at order $`p^6`$ too. On the experimental side, the old radiative correction calculations used in LR  are updated in Radiative  where a clean procedure with generalized form-factors is proposed. The experimental data are mostly analyzed using a linear form factor $`f_+(t)`$. The recent precise CPLEAR data CPLEAR  allow to test this assumption BT2 . Using a linear fit to their data and neglecting systematic errors we get a normalized $`f_+(0)=1`$ and $`\lambda _+=0.0245\pm 0.0006`$. Allowing for curvature we obtain a sizable curvature, $`f_+(0)=1.008\pm 0.009`$ and $`\lambda _+=0.0181\pm 0.0068`$. The fitted curvature is compatible with zero, the central value is precisely at the ChPT prediction given below. In order to obtain $`|V_{us}|`$ with an error of 1% it is therefore important to include the effect of curvature in the analysis. The central value of $`\lambda _+`$ is outside the errors quoted for the linear fit CPLEAR . After the meeting we obtained similar conclusions for the KEK-PS E246 data KEKPS , see v2 of BT2 .

## 4 $`f_+(t)`$ and $`f_0(t)`$: theory

The $`f_+(t)`$ ChPT calculation for is rather cumbersome. I will use our work, BT2  but an independent calculation exists Post3  and agrees reasonably well, the difference is discussed in BT2 . We write the amplitude as

$`f_+(t)`$ $`=`$ $`1+f_+^{(4)}(t)/F_\pi ^2+f_+^{(6)}(t)/F_\pi ^4\text{with}f_+^{(4)}(t)=tL_9^r/2+\text{loops},`$
$`f_+^{(6)}(t)`$ $`=`$ $`8\left(C_{12}^r+C_{34}^r\right)\left(m_K^2m_\pi ^2\right)^2+tR_{+1}^{K\pi }+t^2(4C_{88}^r+4C_{90}^r)+\text{loops}(L_i^r).`$ (4)
The pion electromagnetic form factor data yield BT1  $`L_9^r=0.00593\pm 0.00043`$ and $`4C_{88}^r+4C_{90}^r=0.00022\pm 0.00002.`$ With this input we fit the CPLEAR data and obtain

$$R_{+1}^{K\pi }=(4.7\pm 0.5)10^5\text{GeV}^2\text{and}\lambda _+=0.0170\pm 0.0015.$$
(5)
The first agrees with the VMD estimate BT2  $`R_{+1}^{K\pi }|_{VMD}410^5\text{GeV}^2`$. The latter comes from ChPT as

$$\lambda _+=0.0283(p^4)+0.0011(\text{loops }p^6)0.0124(C_i^r).$$
(6)
Our main conclusion follows from rewriting the full $`p^6`$ result for $`f_0(t)`$ as

$`f_0(t)`$ $`=`$ $`1(8/F_\pi ^4)\left(C_{12}^r+C_{34}^r\right)\left(m_K^2m_\pi ^2\right)^2+(8t/F_\pi ^4)\left(2C_{12}^r+C_{34}^r\right)\left(m_K^2+m_\pi ^2\right)`$ (7)
$`+t/\left(m_K^2m_\pi ^2\right)\left(F_K/F_\pi 1\right)(8t^2/F_\pi ^4)C_{12}^r+\overline{\mathrm{\Delta }}(t)+\mathrm{\Delta }(0).`$
$`\overline{\mathrm{\Delta }}(t)`$ and $`\mathrm{\Delta }(0)`$ contain NO $`C_i^r`$ and only depend on the $`L_i^r`$ at order $`p^6`$ thus ALL needed parameters can be determined experimentally.

$`\mathrm{\Delta }(0)=0.0080\pm 0.0057[\text{loops}]\pm 0.0028[L_i^r],`$ is known and an expression for $`\mathrm{\Delta }(t)`$ can be found in Ref. BT2  The errors are an estimate of higher orders and using fits of the $`L_i^r`$ using different assumptions. The $`p^6`$ estimate of LR  corresponds to $`(8/F_\pi ^4)\left(C_{12}^r+C_{34}^r\right)(m_K^2m_\pi ^2)^20.016\pm 0.008`$.

## 5 Conclusions

I discussed the $`K_\mathrm{}3`$ form factors in ChPT to order $`p^6`$. The main conclusions are that the curvatures for $`f_+(t)`$ and $`f_0(t)`$ can be predicted from the data on pion electromagnetic BT1  and scalar BD  form-factors, the curvature in $`f_+(t)`$ and $`f_0(t)`$ should be taken into account in new precision experiments but from the slope and the curvature we can determine experimentally the needed parameters to calculate $`f_+(0)`$. A precision of better than one percent seems feasible for $`|V_{us}|`$.

This work has been funded in part by the Swedish Research Council and the European Union RTN network, Contract No. HPRN-CT-2002-00311 (EURIDICE)


text:
# Geometry and Conductance of Al Wires Suspended between Semi-Infinite Crystalline Electrodes

## Abstract

We present a first-principles study of a coherent relationship between the optimized geometry and conductance of a three-aluminum-atom wire during its elongation process. Our simulation employs the most definite model including semi-infinite crystalline electrodes using the overbridging boundary-matching method \[Phys. Rev. B 67, 195315 (2003)\] extended to incorporate nonlocal pseudopotentials. The results that the conductance of the wire is $``$ 1 G0 and the conductance trace as a function of electrode spacing shows a convex downward curve before breaking are in agreement with experimental data.

preprint: PRESAT-7803

Atomic-scale wires, which can be produced by scanning tunneling microscopy and break-junction techniques, are the objects of considerable attention datta-ruitenbeek , since the properties of such minute systems with low dimensionality can be different from those of bulks; as the size of the structure becomes smaller than the electronic mean free path, the electron transport becomes ballistic and the conductance is dominated by the electronic structure. To date, there have been a great deal of geometrical for revealing the mechanical and electrical behavior of such minute wires. Among the most exciting discoveries was that a single strand of gold atoms suspended between electrodes manifests conductance quantization in steps of G0=2$`e^2/h`$ rubio-takayanagi-kizuka  as electrode spacing is increased, where $`e`$ is the electron charge and $`h`$ is Planck’s constant. Furthermore, other studies showed that an atomic-scale wire made of aluminum atoms exhibits more complicated and even more interesting features; that is, the plateaus of conductance steps are not flat but have a positive slope before the wire breaks scheer-krans-rubio ; sakai . This somewhat counterintuitive result implies that the conductance increases even though the electrodes are pulled apart.

On the theoretical side, there are many first-principles studies for examining electronic structure and transport properties of atomic-scale wires using tight-binding approximation and/or the structureless jellium electrodes lang1 ; lang2 ; kobayashi1 ; kobayashi2 ; okamoto ; tsukamoto ; wan ; brandbyge ; palacios , some of which reported results consistent with experimental data. For more strict theoretical predictions on details of electron conduction properties, it would be mandatory to completely eliminate artifacts and employ a realistic model where electrodes consist of truly semi-infinite crystals. Recently, the overbridging boundary-matching (OBM) formalism fujimotonano ; fujimoto , which is a novel first-principles treatment of electron transport properties for nanoscale junctions sandwiched between bona-fide semi-infinite crystalline electrodes, has been developed.

In this Letter, we implement elaborate first-principles calculations with the incorporation of the OBM formalism to elucidate a close relationship between the relaxed geometrical structure and electronic conductance of the 3-aluminum-atom wire during the elongation process using a model in which the wire is suspended between two semi-infinite Al(001) bulks. To the best of our knowledge, this is the first theoretical investigation of this relationship where the optimized aluminum wire is connected to genuine semi-infinite aluminum crystalline electrodes which are gradually separated. We find that this 3-aluminum-atom wire has a conductance of $``$ 1 G0, and there is only one channel that gives a dominant contribution to the electron conduction at the Fermi level. The conductance trace as a function of electrode spacing decreases in the early stage of the elongation process and shows a minimum; in the case of shorter electrode spacing, a loop current (LC) is induced around the center atom of the bent wire so as to boost electron transmission. With an increase in electrode spacing, the LC weakens and the conductance reduces. When electrode spacing is still further increased, the wire forms into a straight structure with a one-dimensional character, which contributes to the enhancement of electron transmission. As a result, the curve of the conductance trace is convex downward before the wire breaks, and subsequently the conduction channels between the two electrodes disappear.

First, we determine the optimized atomic structure of the 3-aluminum-atom wire suspended between the Al(001) electrodes during its elongation. Our first-principles simulations are based on the real-space finite-difference method rsfd , which enables us to determine both the self-consistent electronic ground state and the optimized geometry with a high degree of accuracy by virtue of the timesaving double-grid technique tsdg  and the direct minimization of the energy functional dmef . The norm-conserving pseudopotentials ncps  of Troullier and Martins tmpp  are used to describe the electron-ion interaction. Exchange-correlation effects are treated by the local-density approximation lda  in the density-functional theory. The nine-point finite-difference case, i.e., the case of $`N_f`$=4 in Eq. (A2) (here and hereafter, the equation number refers to that in Ref. fujimoto ), is adopted for the derivative arising from the kinetic-energy operator in the Kohn-Sham equation. We set a cutoff energy of 25 Ry, which corresponds to a grid spacing of 0.63 a.u., and further take a higher cutoff energy of 224 Ry in the vicinity of the nuclei with the augmentation of double-grid points tsdg .

For optimization, a conventional supercell, which is indicated by the rectangle enclosed with broken lines in Fig. 1, is employed. The size of the supercell is chosen as $`L_x`$=$`L_y`$=15.12 a.u. and $`L_z`$=22.68 a.u. + 2$`d`$. Here, the $`z`$ axis is taken along the wire axis, the $`x`$ and $`y`$ axes are perpendicular to the wire, $`L_x`$, $`L_y`$ and $`L_z`$ are the lengths of the supercell in the $`x`$, $`y`$ and $`z`$ directions, respectively, and $`d`$ is the average value for the projections of interatomic distances between adjacent atoms of the wire onto the $`z`$ component. The wire is connected to the square bases with the side length of $`\frac{a_0}{\sqrt{2}}`$ at both its ends, modeled after the  aluminum strands, and all of these components intervene between the electrodes produced from three atomic layers of the Al(001) surface. The distance between the electrode surface and the basis as well as that between the basis and the edge atom of the wire is chosen to be 0.5$`a_0`$, where $`a_0`$ (=7.56 a.u.) is the lattice constant of the aluminum crystal. The position of the central atom of the wire is optimized on the (110) plane while the other atoms are fixed during the first-principles structural optimization. We find that the wire exhibits a bent structure below electrode spacing $`L_{es}`$ of 26.46 a.u. When elongated up to $`L_{es}`$=27.09 a.u., the wire manifests geometrical transition from a bent wire to a straight one, and finally distorts at $`L_{es}`$=28.35 a.u. comment1

We next explore the electron transmission of the wire at zero-bias limit with the Landauer formula G=tr(TT)G0, where T is a transmission matrix. The global wavefunctions for infinitely extended states continuing from one electrode side to the other are evaluated by the OBM method fujimoto . For the determination of the effective potential required in the OBM formalism, the electronic structure of the transition region which is composed of wire, bases, and a couple of atomic layers of electrodes (see Fig. 1) as well as the electronic structure of the aluminum crystalline electrode is calculated self-consistently employing a standard supercell geometry with a three-dimensional periodic boundary. The wavefunctions in the crystalline electrode region are computed by the generalized eigenvalue equation of Eq. (19), and in order to circumvent a numerical instability caused by evanescent waves, the wavefunction-matching formula Eq. (11) at the interfaces between the electrode and transition regions is described in terms of the ratio matrices $`R`$’s, which are defined by Eqs. (25) and (B9). Furthermore, for calculating the conduction properties of the aluminum wires with $`s`$\- and $`p`$-valence orbitals, we here extend the OBM scheme to the case including the nonlocal parts of the pseudopotentials of Troullier and Martins tmpp ; since the nonlocal parts are represented by block matrices in the Kohn-Sham Hamiltonian when the real-space finite-difference method is applied, the nine-point finite-difference formula for the kinetic-energy operator is adopted so that the off-diagonal elements of the finite-difference formula can entirely cover the block matrices of the nonlocal parts. Consequently, the Hamiltonian is described by a block-tridiagonal matrix composed of 4$`N_{xy}`$-dimensional block matrices, where $`N_{xy}`$ is the number of grid points on the $`x`$-$`y`$ grid plane comment2 .

Figure 2 shows the channel transmissions at the Fermi level as a function of $`L_{es}`$. The eigenchannels are investigated by diagonalizing the Hermitian matrix TT kobayashi2 . There are certain features that are common to the undistorted 3-aluminum-atom wire: throughout all the ranges of electrode spacing below $`L_{es}`$=28.35 a.u., the first channel, which has the characters of $`s`$ and $`p_z`$ for a straight wire, is widely open, while the second and third ones, which have the characters of $`p_x`$ and $`p_y`$ and are degenerate for the straight wire, have small transmission $``$ $`1.0\times 10^3`$, and the resultant electronic conductance is $``$ 1 G0. Moreover, one can see that the conductance trace as a function of electrode spacing exhibits a convex downward curve having a minimum at $`L_{es}`$=26.46 a.u. Our result regarding the conductance trace accords with the available experimental data scheer-krans-rubio ; sakai , and in particular, is in excellent agreement with the data recently reported by Mizobata et al. sakai .

To gain an insight into this unique conductance trace, we plot in Fig. 3 the electron-density distributions and current densities of the first channel at the Fermi level for several electrode spacings. The states incident from the left electrode are depicted. In the case of a short electrode spacing such as $`L_{es}`$=24.57 a.u., the channel transmission is $``$ 1, and the first channel current forms an LC rotating around the center atom. Upon elongating electrode spacing, one finds that the LC weakens and consequently, the transmission exhibits a minimum. The electron transport through the wire is closely associated with the behavior of the LC: in fact, as the center atom is displaced upward from the most stable position depicted in Fig. 3, the LC becomes enhanced to increase the transmission. When the wire is stretched up to electrode spacing of 27.09 a.u., where it forms into a straight wire from a bent wire, the transmission recovers $``$ 0.9. At this stage, the LC completely disappears and the current is axially symmetric around the wire axis. In the straight wire, the one-dimensional character strengthens, therefore the first channel transmission at the Fermi level is larger than that of the bent wire of $`L_{es}`$ 25.83 a.u. After the wire’s distortion, there are no conduction channels between the two electrodes.

One of the previous theoretical studies on the conductance of the 3-aluminum-atom wire employing structureless jellium electrodes observed a similar current density kobayashi2 . However, this simulation does not have sufficient accuracy to discuss strictly quantitative properties of electron transport because it was implemented using an inferior pseudopotential made only of local parts and a very restricted model where the bond lengths of the wire are fixed to be 5.4 a.u. For a definite interpretation of the relationship between geometry and conductance of the wire, a highly transferable pseudopotential having nonlocal parts and structural optimization are indispensable. Our result that the conductance is $``$ 1 G0 is in accordance with that of the ab initio tight-binding calculation incorporating the square bases made of aluminum atoms palacios . Compared with the other studies using nonlocal pseudopotentials and structureless jellium electrodes without any atomic layers or bases, most of the channel transmissions and conductance in the present study are smaller than those of the other studies: Lang lang1  and Kobayashi et al. kobayashi1  reported that the conductance of the wire is $``$ 1.5 G0, and moreover, Kobayashi et al. reported the transmissions of the second and third channels are $``$ 0.3. In these calculations adopting jellium electrodes without any atomic layers, the channel transmissions and electronic conductance sensitively depend on direct effects from the jellium electrodes, and therefore the difference between our results and theirs may be mainly attributed to the absence of atomic layers or bases.

In summary, we have studied a close relationship between the optimized geometry and electronic conductance of the 3-aluminum-atom wire suspended between genuine semi-infinite aluminum crystalline electrodes. Our findings are that (1) the 3-aluminum-atom wire sandwiched between two Al(001) electrodes exhibits an electronic conductance of $``$ 1 G0, (2) there is one conduction channel which widely opens at the Fermi level, and (3) the conductance trace as a function of electrode spacing behaves so as to draw a convex downward curve before the wire breaks. In the bent wire with a short electrode spacing of $`L_{es}`$26.46 a.u., the electron transport property of the wire depends on the strength of the LC. On the other hand, in the case of a long electrode spacing of $`L_{es}`$27.09 a.u., the one-dimensional character of the straight wire contributes towards the electron conduction to give rise to a conductance of $``$ 1 G0. This result of the conductance trace is in accordance with experimental data. Work in progress is the study on the electron conduction of the aluminum wire attached to the electrodes having other crystal faces such as the (111) and (110) planes.

This research was partially supported by the Ministry of Education, Culture, Sports, Science and Technology, Grant-in-Aid for Young Scientists (B), 14750022, 2002. The numerical calculation was carried out by the computer facilities at the Institute for Solid State Physics at the University of Tokyo, and Okazaki National Institute.