text:
# Techniques for effective vocabulary selection

## 1 Introduction

The size and performance of a language model or speech recognition system are often strongly influenced by the size of its vocabulary. Ideally, the vocabulary is small, allowing us to build compact language models, and it is matched to the target domain so that as many as possible of the domain-specific words are known to the recognition system. Compact language models generate compact word graphs that are efficient to search and domain-matched vocabularies result in fewer out-of-vocabulary (OOV) words and consequently fewer recognition errors. In a study of the effect of OOV words on the Word Error Rate (WER) of a recognition system, Rosenfeld  arrives at a figure of about 1.2 WER points per OOV word in a typical large vocabulary task.

While a large and comprehensive vocabulary may be desirable from the point of view of lexical coverage, we often settle for smaller and more tractable vocabularies. Not only are large vocabulary language models themselves very large, but for speech recognition systems, there is also the additional cost and effort involved in determining accurate pronunciations for every vocabulary entry. Even with the help of tools to generate pronunciations and check consistency of entries, this is a difficult task .

Furthermore, there is also the attendant problem of increased acoustic confusability for speech recognition systems when the vocabulary is large . For applications requiring a finite vocabulary, picking the right words for the vocabulary is especially important for achieving satisfactory performance. Usually, a number of text corpora from various domains and time periods are available on which to train. The target domain is known, and the amount of data available in the target domain is far less than in any of the training corpora. Clearly, restricting the vocabulary to just the words that are observable in a meager amount of available domain data would be disastrous. On the other hand, including the union of the vocabularies of all the training corpora would be intractable. What we want in this situation is to assume that the vocabulary of the target domain is somehow related to the vocabularies of each training corpus, and subsequently infer the target vocabulary from the individual training corpus vocabularies, considering the observable portion of the domain text to be a sample.

Even though vocabulary selection is an important issue and the problem appears to be simple, little work exists on this topic to date. The most common approaches seem to be ad hoc in nature, typically including words from each corpus that exceed some threshold frequency. This threshold depends on intuitions about the relevance of the corpus to the target domain . In Rosenfeld’s 1995 work  on optimizing vocabularies, attention was mainly directed at determining the effect on the OOV rate of corpus recency, size and origin. While it was found that all three factors strongly affected the OOV rate, no specific guidelines were proposed as to how to combine the vocabularies from these different corpora to choose the target vocabulary. Indeed, Rosenfeld remarks that an ad hoc approach that discounted words by 1% for every week of age of the corpus reduced the OOV rate only very slightly for vocabulary sizes in the range of 20,000 to 50,000 words.

The paucity of work on this important topic can partly be attributed to the general observation due to Zipf  that with even a moderate sized vocabulary chosen wisely, one can hope to get significant lexical coverage. Yet it is desirable from the point of view of scalability, extensibility and generality to study principled methods to address this problem. In this paper, we propose three such principled methods. The goal is to select a single vocabulary from many corpora of varying origins, sizes and recencies such that the vocabulary is optimized for both size and the OOV rate in the target domain. Section 2 defines the problem. Section 3 describes the proposed techniques, and Section 4 presents the results.

## 2 Problem Description

The vocabulary selection problem can be briefly summarized as follows. We wish to estimate the true vocabulary counts of a partially visible corpus of in-domain text (which we call the held-out set) when a number of other fully visible corpora, possibly from different domains, are available on which to train. There is an implicit assumption that the held-out text is related to the training text and the learning task amounts to inferring this relation. The reason for learning the in-domain counts $`x_i`$ of words $`w_i`$ is so that the words may be ranked in order of priority, enabling us to plot a curve relating a given vocabulary size to its OOV rate on the held-out corpus. Therefore, it is actually sufficient to learn some monotonic function $`f\left(x_i\right)`$ in place of the actual $`x_i`$. We may assume that the counts are normalized by document length so that the amount of available data for a particular corpus is itself irrelevant to the issue at hand.

Table 1 illustrates the problem; $`n_{i,j}`$ are the visible counts from each of the documents $`j`$, for the word $`w_i`$, and the $`x_i`$ are the incomplete counts for words $`w_i`$ in the partially observable domain text.

Let $`x_i`$ be some function $`\mathrm{\Phi }_i`$ of the known counts $`n_{i,j}`$ for $`1<=j<=m`$ for each of the $`m`$ corpora. Then, the problem can be restated as one of learning the $`\mathrm{\Phi }_i`$ from a set of examples where

$$x_i=\mathrm{\Phi }_i(n_{i,1},\mathrm{},n_{i,m})$$

In the following section, we summarize three techniques for learning the $`\mathrm{\Phi }_i`$ that optimize the vocabulary for the domain from which the held-out data was drawn.

## 3 Method

For simplicity, let the $`\mathrm{\Phi }_i`$ be linear functions of the $`n_{i,j}`$ and that they are independent of the particular words, $`w_i`$. That is, $`\mathrm{\Phi }=\mathrm{\Phi }_i=\mathrm{\Phi }_j,i,j`$. Then, we can write

$$\mathrm{\Phi }(n_{i,1},\mathrm{},n_{i,m})=\underset{j}{}\lambda _jn_{ij}$$
(1)
The problem transforms into one of learning the $`\lambda _j`$. We now outline three methods to do this. The first is based on maximum likelihood (ML) count estimation, the second and third are based on document similarity measures. We evaluate each of these three methods against a fourth baseline method that simply assigns identical values to all the $`\lambda _j`$.

### 3.1 Maximum likelihood count estimation

In ML count estimation, we simply interpret the normalized counts $`n_{ij}`$ as probability estimates of $`w_i`$ given corpus $`j`$ and the $`\lambda _j`$ as mixture coefficients for a linear interpolation. We try to choose the $`\lambda _j`$ that maximize the probability of the in-domain corpus. Formally, let $`\mathrm{P}\left(w_i|j\right)=n_{i,j}`$. Our goal is to find

$$\widehat{\lambda _1},\mathrm{},\widehat{\lambda _m}=\underset{\lambda _1,\mathrm{},\lambda _m}{argmax}\underset{i=1}{\overset{\left|V\right|}{}}\left(\underset{j}{}\lambda _j\mathrm{P}\left(w_i|j\right)\right)^{C\left(w_i\right)}$$
(2)
where $`C\left(w_i\right)`$ is the count of $`w_i`$ in the partially observed held-out corpus and $`V`$ is the set of words in the vocabulary. The $`\lambda _j`$ can subsequently be estimated via the EM algorithm  and used to calculate the interpolated normalized counts. The procedure shown in Figure 1, for instance, is effective in rapidly computing the values of the $`\lambda _j`$.

### 3.2 Document-similarity-based count estimation

The document-similarity-based count estimation method calculates interpolation weights from similarity measures between the held-out corpus and each of the training corpora. This similarity measure can presumably be calculated using any of a number of methods ranging from a simple Euclidean distance metric to a more sophisticated divergence measure between the observable probability distributions, such as Kullback-Liebler (KL)  or a symmetric variant .

The Euclidean distance metric is calculated as follows. Suppose we represent each document by a vector of its normalized word counts. Then, the Euclidean distance between two corpora $`C_a`$ and $`C_b`$, $`\mathrm{\Delta }(C_a,C_b)`$, is given by

$$\mathrm{\Delta }(C_a,C_b)=\sqrt{\underset{i=1}{\overset{\left|V\right|}{}}\left(n_{a,i}n_{b,i}\right)^2}$$
(7)
where $`n_{a,i}`$ and $`n_{b,i}`$ are the normalized counts of $`w_i`$ in the corpora $`a`$ and $`b`$, respectively.

Likewise, the KL-divergence, which we again denote as $`\mathrm{\Delta }(C_a,C_b)`$ for the sake of uniformity, is given by

$$\mathrm{\Delta }(C_a,C_b)=\underset{i=1}{\overset{\left|V\right|}{}}\mathrm{P}(a,i)\mathrm{log}_2\left(\frac{\mathrm{P}(a,i)}{\mathrm{P}(b,i)}\right)$$
(8)
where we interpret the normalized word counts as probabilities.

In each of the above distance calculation schemes, let $`D_j`$ be the distance of the $`j`$th corpus from the held-out domain text. Then, since the relevance of a corpus to the domain is inversely related to its distance from the domain, we define

$$\lambda _j=\frac{1/D_j}{_k1/D_k}$$

### 3.3 Data sources and implementation

The experimental setup consisted of learning the optimal vocabulary to model the language of the English broadcast news. A small amount of hand-corrected closed captioned data, amounting to just under 3 hours (about 25,000 words), drawn from six half-hour broadcast news segments from January 2001, was used as the partially visible held-out data to estimate the two mixture weights $`\lambda _1`$ and $`\lambda _2`$. This held-out data is part of the corpus released by the Linguistic Data Consortium (LDC) for the National Institute of Standards and Technology (NIST) sponsored English topic detection and tracking (TDT4) task.

The training corpora were deliberately chosen to be as different from each other in character as possible. The first corpus consisted of about 18.5 million words of English newswire data covering the period July 1994 through July 1995, and was distributed by the LDC for the NIST-sponsored Hub3 1995 continuous speech recognition task. It contained text from The NY Times News Service, LA Times, Washington Post News Service, Wall Street Journal and Reuters North American Business News. The second training corpus consisted of a closer match to the target domain and came from segments of the TDT4 dataset released by the LDC. This consisted of about 2.5 million words of closed captioned transcripts from the period November through December 2000.

Unigram counts for the training and held-out corpora were generated using language modeling tools from the SRILM  using Witten-Bell  smoothing. Estimation of the $`\lambda _j`$ was performed on five of the six held-out segments which we collectively refer to as the development corpus, and OOV rates were measured on the remaining segment, which we refer to as the test corpus. This procedure was repeated six times, one for each possible split of the held-out data. The results we present are averaged numbers obtained from the six splits. Where applicable, we use the subscripts “hub3” and “tdt4” to refer to parameters specific to the above corpora.

## 4 Results and Discussion

We examine the results of our experiments to evaluate the various methods. Figure 2 shows a plot of the OOV rate against increasing vocabulary size from 1 word to 90,000 words. This figure, which is plotted in the logarithmic scale, is only meant to show the general shape of the individual plots and for drawing some broad generalizations. For instance, we see confirmation of the common observation that the OOV rate of a given vocabulary on a corpus is logarithmically related to the vocabulary size. Furthermore, it is also evident that for small vocabularies there exist obvious differences in the performances of the various vocabulary selection methods. But for large vocabularies, this difference is seen to diminish. Indeed, for vocabulary sizes in excess of about 60,000 words, the four plots practically merge into a single line showing that at around that threshold and beyond, we capture practically all the words that are likely to be used in the domain under consideration, regardless of the specific method used to choose the vocabulary.

For a finer-grained comparison of the individual techniques, we restrict our attention to the rectangular sub-region in Figure 2, which is depicted in a separate plot in Figure 3. It shows the performance of the four systems for a vocabulary range of 1,000 to 2,000 words. The trend of the curves in this graph, which continues up to a vocabulary size of around 40,000 words, clearly shows that the ML method outperforms all the other three methods by over 1% absolute. It is also clear that the method based on KL-divergence is the poorest of all, performing worse than even the uniform baseline. The Euclidean-distance-based method performs almost identically as the uniform baseline (and thus the plot for the latter, being almost hidden behind that of the former, is barely visible).

In hindsight, the relatively good performance of the maximum-likelihood-based method is not very surprising because it is the only method that does not look beyond the development corpus vocabulary to compute its objective function. Both the KL-divergence-based method and the Euclidean-distance-based method sum quantities over the entire vocabulary and are therefore affected by the values held by individual words that were not seen in the development corpus. This problem is especially acute because the actual vocabulary of the partially visible development corpus is typically tiny compared to the vocabularies of the training corpora. The KL-divergence-based method is affected most by this situation. Because KL-divergence involves calculation of log-probabilities, the method is extremely sensitive to the amount of probability mass devoted to unseen vocabulary items and consequently to the particular form of smoothing employed. Since a significant number of words in the vocabulary are typically unseen in the development corpus, these end up with very low unigram probabilities. Thus, in summing over the entire vocabulary, large negative numbers come into play which overshadow any significant contribution to the total divergence by the unigrams observed in the development corpus. We suspect therefore that we must not attach much significance to the final quantity computed by this method unless the size of the development corpus itself is substantial.

The Euclidean method is also likewise affected, but to a lesser degree and slightly differently. The computed distances tend to be dominated by words that are absent in the development corpus rather than by words that are present in it. Since the absent words form the bulk of the vocabulary, the distances computed between the various corpora and the development text, and consequently the $`\lambda _j`$ will all roughly be the same, as evidenced by the figures in Table 2.

## 5 Conclusions

We have outlined three general techniques to select an optimal vocabulary for domain-specific speech and language modeling tasks. The techniques are scalable to arbitrarily large-sized corpora and extensible to any number of corpora. Whenever reasonable amounts of training data and reliable unigram count estimates are available, we believe that the maximum-likelihood-based method we have described is a robust way to select a domain’s vocabulary especially when its size is expected to be under a certain threshold. This threshold can be expected to vary between domains and it is possible that when it is high, the choice of any particular strategy over another does not matter. However, we believe that always following a principled strategy to select the vocabulary offers the safest path.

We plan to continue to refine and evaluate the techniques presented in this paper and apply them for vocabulary selection in the English broadcast news recognition task of the NIST 2003 Hub4 evaluation.

## 6 Acknowledgments

This material is based upon work supported by the Defense Advanced Research Projects Agency under Contract No. MDA972-02-C-0038. The authors also thank Andreas Stolcke and Kristin Precoda for useful discussions regarding its content and Judith Lee for careful proof-reading of the manuscript.


format: text

text:
# Proposal for a simple quantum error correction test gate in linear optics

## I Introduction

Whilst high hopes are held for the eventual demonstration of large scale quantum processing, present experimental attention remains fixed on few qubit demonstrations. Implementation of few qubit algorithms in the competing platforms provides valuable insight into the important physics and technical issues of different architectures.

A relatively new contender for scalable quantum computation is the scheme due to Knill et al KLM  based on single photon, dual rail qubits; linear optical networks; and photon resolving measurement and feedforward. Simple gates based on this scheme have been demonstrated pit02  and more general gates are planned ral01 ; pit01 ; ral02 . It is thus timely to consider what small scale circuits might be possible with the currently available technology.

A key enabling quantum circuit is error correction shor95 ; ste96 . Even medium scale quantum processing is expected to be impossible without error correction. Here we describe a simple example of error correction which would suit demonstration with a linear optical circuit. In its simplest form the experiment would require only three coincident photons.

## II The Circuit

The quantum circuit we wish to consider is a simplification of the standard bit flip correcting code nie00  and is shown in Fig.1. The qubit is encoded with the following relationship between logical, $`|_L`$, and physical, $`|`$, qubits:

$`|0_L`$ $`=`$ $`|0|0`$
$`|1_L`$ $`=`$ $`|1|1`$ (1)
We suppose that one of the physical qubits, say the second, suffers decoherence which produces random bit-flips. As a result of this decoherence an arbitrary initial qubit:

$`\alpha |0_L+\beta |1_L=\alpha |0|0+\beta |1|1`$ (2)
evolves into the mixed state:

$`\rho `$ $`=`$ $`(1P)(\alpha |0|0+\beta |1|1)(\alpha ^{}0|0|+\beta ^{}1|1|)+`$ (3)
$`P(\alpha |0|1+\beta |1|0)(\alpha ^{}0|1|+\beta ^{}1|0|)`$
where $`P`$ is the probability that the bit flip occurred. The aim of the circuit which follows is to return the qubit to its original logical value. To achieve this requires an ancilla photon prepared in the physical zero state. We can write the combined state of the qubits and ancilla after decoherence as

$`\rho `$ $`=`$ $`(1P)|0_a(\alpha |0|0+\beta |1|1)(\alpha ^{}0|0|+\beta ^{}1|1|)0|_a+`$ (4)
$`P|0_a(\alpha |0|1+\beta |1|0)(\alpha ^{}0|1|+\beta ^{}1|0|)0|_a`$
with the subscript “a” labelling the ancilla state. A controlled not (CNOT) gate is applied with the ancilla as target and the first qubit as control. This transforms the state to

$`\rho `$ $`=`$ $`(1P)(\alpha |0|0|0_a+\beta |1|1|1_a)(\alpha ^{}0|0|0|_a+\beta ^{}1|1|1|_a)+`$ (5)
$`P(\alpha |0|1|0_a+\beta |1|0|1_a)(\alpha ^{}0|1|0|_a+\beta ^{}1|0|1|_a)`$
A second CNOT is then applied with the ancilla as target but now the second qubit acts as control. The state becomes

$`\rho `$ $`=`$ $`(1P)|0_a(\alpha |0|0+\beta |1|1)(\alpha ^{}0|0|+\beta ^{}1|1|)0|_a+`$ (6)
$`P|1_a(\alpha |0|1+\beta |1|0)(\alpha ^{}0|1|+\beta ^{}1|0|)1|_a`$
Finally we detect the ancilla state. If we find the ancilla in the zero state then the logical qubit is projected onto its original state and no correction is neccessary. On the other hand if the ancilla is found in the one state then the projected state is

$`\alpha |0|1+\beta |1|0`$ (7)
We know an error has occurred which can be corrected by flipping the value of the second physical qubit and thus returning the logical qubit to its initial value.

The circuit can also be understood in the language of stabilizer codes got96 . The code space is the $`+1`$ eigenstates of $`ZZ`$, whilst the error space is the $`1`$ eigenstates of $`ZZ`$. Here $`Z`$ is the Pauli sigma z operator and $`ZZ`$ indicates the tensor product of a sigma z measurement on the first qubit with a sigma z measurement on the second qubit. The two CNOT’s achieve precisely this measurement with the ancilla equals zero result indicating the $`+1`$ eigenstate, and hence the system being in the code space, while the ancilla equals one result indicates the $`1`$ eigenstate and hence the system is in the error space and needs to be corrected.

The usefulness of this circuit is obviously limited by the very specific nature of the errors corrected, ie only bit flips on one of the physical qubits. Never-the-less it exhibits the same basic structure as more versatile codes nie00  whilst limiting the complexity of the required circuit. In the next section we will discuss an optical implementation of this code which appears tractable to current experimentation.

## III The Optical Implementation

We consider single photon qubits with the polarization degree of freedom determining their logical values. We define the physical qubit value “zero” as being a single horizontally polarized photon, $`|0|H`$, and a physical qubit value of “one” as being a single vertically polarized photon, $`|1|V`$. Knill et al KLM  showed that non-deterministic CNOT gates could be constructed from linear optics with success rates of one in sixteen, using two additional ancilla photons. A CNOT would be required to produce the encoded logical qubits. Thus an implementation of our error correction circuit would require three CNOTs and hence require the simultaneous production of nine single photon states. The success rate would be about one in four thousand. Although not beyond the realm of medium term possibility, such an experiment is currently not feasible. However in the following we will discuss an in principle demonstration utilizing the coincidence basis with much lower technical requirements.

The proposed set up is shown schematically in Fig.2. The encoded state is produced directly using type one down conversion through a pair of $`\chi 2`$ crystals with orthogonally oriented optical axes whi99 . Pairs of photons originating from one crystal will be horizontally polarized whilst those originating from the other crystal will be vertically polarized. By changing the polarization of the pump beam the proportion of horizontal to vertical pairs can be continuously varied. In the far field, for sufficiently thin crystals, spatial information on the origin of the pairs is erased and the resulting entangled output state is approximately

$`|vac_1|vac_2+\chi (\alpha |H_1|H_2+\beta |V_1|V_2)`$ (8)
where $`|vac_i`$ is the vacuum state, $`1`$ and $`2`$, label the two beams, $`\alpha =\mathrm{cos}\theta `$ and $`\beta =\sqrt{1\alpha ^2}`$ and $`\theta `$ is the orientation of the pump beam polarization away from vertical. Coincidence detection will pick out only the doubly occupied parts of the state, so the effective input state is

$`\alpha |H_1|H_2+\beta |V_1|V_2`$ (9)
which is logically identical to Eq.2.

Controlled decoherence can be introduced onto the second qubit by passing it through a birefringent crystal oriented at $`45`$ degrees to horizontal jam01 . The effect of the crystal is to pull apart in time the two polarization modes in the diagonal/anti-diagonal basis. When the time difference becomes an appreciable fraction of the coherence length decoherence occurs. In the horizontal/vertical basis the resut is random bit flips and the state produced can be written in the form of Eq.3. The degree of decoherence, $`P`$, is a function of the crystal length.

To detect the errors we must introduce an additional ancilla qubit for readout. This can be supplied by a second down-converter producing just horizontal pairs. A CNOT gate which works in the coincidence basis and does not require additional ancilla modes can be implemented using just linear optics ral02 ; hof02 . The schematic of Fig.2 shows the optical network needed to implement the two CNOT’s in the quantum circuit (Fig.1). The polarization qubits are decomposed into separate spatial modes using polarizing beamsplitters then the modes from different qubits are mixed on beamsplitters. Many optical paths are possible through the network but only a few result in photonic qubits at all three outputs, as determined by coincident detection of photons. In such cases quantum interference due to indistinguishability of photons ensures the required transformations are implemented. The success rate is one in eight-one.

Correction of the decohered qubit could be implemented, if required by the result of the ancilla detection, using a fast Pockel cell pit02b .

For most runs of the experiment insufficient photons will be detected and we have a null result. However, on those occasions when a photon is detected at the ancilla trigger and at the ancilla output, and photons are detected at both qubit outputs, then to an excellent approximation the quantum circuit of Fig.1 will have been implemented. For such event we would expect the logical value of the qubits should ideally be the same as that prepared, in spite of the presence of the decohering element.

Four photon coincidences of a few per minute have been achieved with down conversion bou02 . This will be further reduced by a factor of one in eighty-one for this proposal. On the other hand working in the coincidence basis means, at least in principle, that this source efficiency and the efficiency of the detectors does not effect the fidelity of the accepted events.

## IV Conclusion

We have described a simple error correction code, and proposed an in principle test of its operation using current four-photon technology. Although both the code and its implementation are major simplifications over what would be required in a scalable architecture, many of the basic principles are common. We thus suggest that pursuit of experiments like the one proposed here will reveal much about the important physical and technical issues to be faced for truly scalable architectures.

###### Acknowledgements.

We thank Charlene Ahn and Gerard Milburn for motivating discussions. This work was supported by the Australian Research Council and ARDA.


text:
# Quantum tunneling in a three dimensional network of exchange coupled single-molecule magnets

## Abstract

A Mn4 single-molecule magnet (SMM) is used to show that quantum tunneling of magnetization (QTM) is not suppressed by moderate three dimensional exchange coupling between molecules. Instead, it leads to an exchange bias of the quantum resonances which allows precise measurements of the effective exchange coupling that is mainly due to weak intermolecular hydrogen bounds. The magnetization versus applied field was recorded on single crystals of \[Mn4\]2 using an array of micro-SQUIDs. The step fine structure was studied via minor hysteresis loops.

Single-molecule magnets (SMM), such as Mn12, Mn4 and Fe8 Christou00 ; Sessoli93b ; Sessoli93 ; Aubin96 ; Boskovic02 , have become model systems to study quantum tunneling of magnetization (QTM) Friedman96 ; Thomas96 ; Sangregorio97 ; Aubin98 ; Aubin98b ; WW\_Science99 . These molecules comprise several magnetic ions, with their spins coupled by strong exchange interactions to give a large effective spin. The molecules are regularly assembled in large crystals where often all the molecules have the same orientation. Hence, macroscopic measurements can give direct access to single molecule properties. Many non-magnetic atoms surround the magnetic core of each molecule. Exchange interactions between molecules are therefore relatively weak and have been neglected in most studies.

Recently, the study of a dimerized SMM \[Mn4\]2 showed that intermolecular exchange interactions are not negligible WW\_Nature02 . This compound belongs to the \[Mn4O3Cl4(O2CR)3(py)3\]2 family, with R = CH2CH3 and it will be named in the following as compound 1. The spins of the two Mn4 molecules are coupled antiferromagentically. Each molecule acts as a bias on its neighbor, the quantum tunneling resonances thus being shifted with respect to the isolated SMM. In this letter we show that even in three-dimensional networks of exchange coupled SMMs, ordering effects do not quench tunneling.

We selected a dimerized SMM \[Mn4\]2, called compound 2. The molecule belongs to the same family as compound 1, however R = CH3. Because this substituent has a smaller volume than R = CH2CH3, molecules are packed closer together. This leads to stronger interdimer interactions, which are negligible in compound 1. The preparation, X-ray structure and detailed physical characterization have been reported elsewhere Hendrickson92 ; Aliaga\_Alcade . Compounds 1 and 2 crystallize in the hexagonal space group $`R3`$(bar) with two Mn4 molecules per unit cell lying head-to-head on a crystallographic S6 symmetry axis (Fig. 1). Each monomer Mn4 has a ground state spin S = 9/2. The Mn-Mn distances and the Mn-O-Mn angles are similar and the uniaxial anisotropy constant is expected to be the same for the two dimer systems. These dimers are held together via six C$``$H$`\mathrm{}`$Cl hydrogen bonds between the pyridine (py) rings on one molecule and the Cl ions on the other and one Cl$`\mathrm{}`$Cl Van der Waals interaction (Fig. 1a). These interactions lead to an antiferromagnetic superexchange interaction between the two Mn4 units of a dimer WW\_Nature02 .

Owing to the S6 symmetry of \[Mn4\]2, each Mn4 is close to three neighboring Mn4 molecules of the three neighboring \[Mn4\]2 (Fig. 1b). There are hydrogen bonds between the pyridine (py) rings of the molecules and the O ions of the other three neighboring molecules. The C$``$H$`\mathrm{}`$O hydrogen bonds between \[Mn4\]2 dimers have C$`\mathrm{}`$O distances and C$``$H$`\mathrm{}`$O angles of 4.05 Å and 124.85, respectively. The interactions between the dimers are expected to be antiferromagnetic and weaker than the interdimer interactions. The two different antiferromagnetic couplings, the stronger one inside the dimer and the weaker one between the dimers, make this system an interesting candidate for studying the QTM in a three dimensional network of exchange coupled SMMs.

The magnetization versus applied field was recorded on single crystals of \[Mn4\]2 using an array of micro-SQUIDs WW\_ACP\_01 . Figs. 2a and 2b show typical hysteresis loops of magnetization versus applied field for different field sweep rates. The field is applied along the easy axis of magnetization of a single crystal of about 20 $`\mu `$m. These loops display step-like features separated by plateaus. The hysteresis loops of the two crystals are similar. However, compound 2 shows a fine structure that is absent in the hysteresis loops of compound 1. We will show in the following that the main feature of the hysteresis loops can be explained by the QTM of one Mn4 molecule, coupled by a superexchange interaction $`J`$ to the other unit of the \[Mn4\]2 dimer. We discuss first compound 1 because the coupling with neighboring dimers can be neglected WW\_Nature02 . Then, we show that the fine structure observed for compound 2 is induced by a superexchange interaction $`J^{}`$ between neighboring dimers (Fig. 1b).

The simplest Hamiltonian describing the spin system of an isolated SMM is:

$$=DS_z^2+_{\mathrm{trans}}+g\mu _\mathrm{B}\mu _0\stackrel{}{S}\stackrel{}{H}$$
(1)
$`S_x`$, $`S_y`$, and $`S_z`$ are the components of the spin operator; $`D`$ is the anisotropy constant defining an Ising type of anisotropy; $`_{\mathrm{trans}}`$, containing $`S_x`$ or $`S_y`$ spin operators, gives the transverse anisotropy which is small compared to $`DS_z^2`$ in SMMs; and the last term describes the Zeeman energy associated with an effective field $`\stackrel{}{H}`$. For one isolated spin the effective field is the applied field. This Hamiltonian has an energy level spectrum with $`(2S+1)`$ values which, to a first approximation, can be labelled by the quantum numbers $`M=S,(S1),\mathrm{},S`$ taking the $`z`$-axis as the quantization axis. The energy spectrum can be obtained by using standard diagonalization techniques. At $`\stackrel{}{H}=0`$, the levels $`M=\pm S`$ have the lowest energy. When a positive field $`H_z`$ is applied, the levels with $`M>0`$ decrease in energy, while those with $`M<0`$ increase. Therefore, energy levels of positive and negative quantum numbers cross at certain values of $`H_z`$ given by $`\mu _0H_znD/g\mu _\mathrm{B}`$, where $`n=0,1,2,3,\mathrm{}`$. When the spin Hamiltonian contains transverse terms ($`_{\mathrm{trans}}`$), the level crossings can be avoided level crossings. The spin $`S`$ is in resonance between two states when the local longitudinal field is close to an avoided level crossing. The energy gap, the so-called tunnel splitting $`\mathrm{\Delta }`$, can be tuned by a transverse field (a field applied perpendicular to the $`z`$ direction) via the $`S_xH_x`$ and $`S_yH_y`$ Zeeman terms. The effect of these avoided level crossings leads to well defined steps in hysteresis loop measurements.

The main point to note is that the giant spin Hamiltonian predicts always the first level crossing at zero field, corresponding to the QTM of a SMM between $`M=\pm S`$ states. Thus, for compound 1 (see Fig. 2a), the single-spin Hamiltonian is not sufficient to explain the first resonance shifted to negative fields and the absence of the quantum tunnelling at zero field, in contrast to other SMMs.

In order to explain the observed features in Fig. 2a, one has to take into account the superexchange coupling $`J`$ between pairs of Mn4 units. A Hamiltonian for the two-coupled molecules can be written and the energy states of the \[Mn4\]2 can be calculated by exact diagonalization. More details on the dimer Hamiltonian and the corresponding Zeeman diagram are reported elsewhere WW\_Nature02 . Here, we propose a phenomenological model that is sufficiently simple to allow inclusion, in a second step, of more exchange couplings. The influence of the exchange coupling of the neighboring molecule is taken into account by an exchange bias field $`H_{bias}`$. The effective field $`H_z`$ acting on the molecule is therefore the sum of the applied field $`H_{app}`$ and the bias field $`H_{bias}`$ :

$$H_z=H_z^{\mathrm{app}}+H_z^{\mathrm{bias}}=H_z^{\mathrm{app}}+\frac{J}{g\mu _\mathrm{B}\mu _0}M_2$$
(2)
where $`M_2`$ is the quantum number of the neighboring molecule and $`J`$ is the associated exchange coupling. In the following we explain the hysteresis loops when the field $`H_z^{\mathrm{app}}`$ is swept from negative to positive values. At low temperature, $`M_2`$ has two possible values $`M_2=\pm S=\pm 9/2`$. We expect therefore resonant QTM for applied fields $`\mu _0H_z^{\mathrm{app}}nD/g\mu _\mathrm{B}\pm M_2J/g\mu _\mathrm{B}`$, where $`n=0,1,2,3,\mathrm{}`$. The two possibilities of $`M_2`$ are represented by two combs in Fig. 2a. The first comb (0,1,2) corresponds to $`M_2=9/2`$ and the second one (0’,1’) to $`M_2=9/2`$. This model describes all observed quantum transitions in Fig. 2a with two fitting parameters $`D/k_\mathrm{B}=0.72`$ K and $`J/k_\mathrm{B}=0.1`$ K. It neglects co-tunneling and other two-body tunnel transitions having a lower probability of occurrence WW\_Nature02 ; WW\_PRL02 .

Compound 2 displays hysteresis loops (Fig. 2b) similar to those of compound 1. However, the total exchange coupling is larger for compound 2. The values of $`D/k_\mathrm{B}=0.75`$ K and $`J/k_\mathrm{B}`$ = 0.15 K were obtained from the field positions of the steps in the hysteresis loops. Another difference between the two compounds is that the hysteresis loops of compound 2 exhibit fine structure that can not be explained by the dimer model described above (Eq. 2). In order to better analyze this fine structure, minor hysteresis loops were measured (Figs. 3 and 4). First the sample is saturated in positive field; all the molecules are in the $`M=+9/2`$ state. Then the field is decreased. The system approaches the first avoided energy level crossing at a field value of $``$0.5 T. A fraction of the dimers switches from +9/2 to $``$9/2, and the total magnetization of the system decreases, generating a step in the hysteresis loop. When the magnetization reaches the second plateau ($``$0.2 T), the field is swept back towards positive saturation; the tunneling from $`M=9/2`$ to 9/2 is favored via the exited state 7/2 ($``$1 T). After this transition the sample reaches positive saturation. The purpose of these minor hysteresis loops is to confirm the fine structure of each transition starting from different initial states.

The tunnel transitions exhibit four equidistant kinks, that we explain by the exchange coupling to the three neighboring dimers. The spin of the three neighboring molecules can be either aligned with the magnetic field or reversed, leading to four different situations: from zero to three reversed neighbors.

The exchange coupling between a molecule and its neighbors acts like a supplementary field bias and shifts further the resonance fields. The total field bias induced by the neighbors and the other Mn4 unit of the dimer can be written:

$$H_{\mathrm{bias}}^{tot}=\frac{1}{g\mu _B\mu _0}\left(JM_2+\underset{i=1}{\overset{3}{}}J^{}M_i^{}\right)$$
(3)
where the first term is the contribution of the intradimer coupling, and $`M_i^{}`$ is the quantum number of the three neighboring dimer molecules $`i`$ (Fig. 1b).

After positive saturation all the molecules are aligned with the field. The first kink in the hysteresis loop corresponds to the QTM of one molecule in the bias field of its non-reversed neighbors. The resonance is shifted towards negative values by the bias field $`H_{\mathrm{bias}}=9/(2g\mu _B\mu _0)(J+3J^{})`$ (see Eq.4). After this first kink, some molecules now have one reversed neighbor. At the second kink it is this newly created population which tunnels generating molecules with two reversed neighbors. The corresponding field shift is $`H_{\mathrm{bias}}=9/(2g\mu _B\mu _0)(J+J^{})`$. The third and the fourth kinks are generated by the QTM of molecules having, respectively, two and three reversed neighbors. The field shift between two consecutive kinks is $``$0.1 T, corresponding to an interdimer interaction $`J^{}`$ 0.015 K.

Minor hysteresis loops were measured for different field sweep rates (Fig. 3) and reversal fields (Fig. 4) in order to probe the step heights of the fine structure: the smaller the sweep rate the higher the resulting kink. This dependence is justified by the Landau Zener model. The main point to note is that heights of two consecutive kinks are correlated. The second kink height is smaller than the first kink height, the third smaller than the second, and so on. This result is in good agreement with our model: in order to have quantum tunneling of molecules with $`n`$ reversed neighbors, the $`n`$ neighbors must have previously reversed.

All the other transitions exhibit the same kind of fine structure, which can be explained by the above model leading to the 8 combs in Fig. 2b, giving for the three fitting parameters $`D/k_\mathrm{B}0.75`$ K, $`J/k_\mathrm{B}0.1`$ K and $`J^{}/k_\mathrm{B}0.015`$ K.

The above results demonstrate that a three dimensional network of exchange coupled SMMs doesn’t suppress QTM. The intermolecular interactions are strong enough to cause a clear field bias, but too weak to transform the spin network into a classical antiferromagnetic material. This three dimensional network of exchange coupled SMMs demonstrate that the QTM can be controlled using exchange interactions, and opens up new perspectives in the use of supramolecular chemistry to modulate the quantum physics of these molecular nanomagnets.


text:
# Structural, electronic, and magnetic properties of a ferromegnetic semiconductor: Co-doped TiO2 rutile

## I INTRODUCTION

Spintronics has emerged expeditiously as a novel technology where both charge and spin are harnessed for new functionalities in devices combining standard microelectronics with spin-dependent effect.Ohno (1998); Ohno et al. (1999); Wolf et al. (2001). For their practical spintronic applications, a ferromagnetic semiconductor is required to be not only easily integerable but also having a high Curie temperature. The development of new materials grows rapidly. Room-temperature ferromagnetism has been observed, for instance, in Mn doped ternary compounds such as CdGeP2Medvedkin et al. (2000), ZnSnAs2Choi et al. (2002), and ZnGeP2Cho et al. (2002). A recent discovery of room-temperature ferromagnetism in Co-doped TiO2 anatase by Matsumoto $`et`$ $`al`$. using combinatorial molecular beam epitaxy (MBE) techniqueMatsumoto et al. (2001) has motivated intensive studies on the structural and physical properties of this material.Chambers et al. (2001); Park et al. (2002a); Kim et al. (2002); Geng and Kim (2002) Transparent to visible light, this material is expected to be of particular importance in spin-based optoelectronic applications such as spin-light-emission-diode and optical switches operating at terahertz frequencyWolf et al. (2001).

It is interesting to note that Matsumoto $`et`$ $`al`$. observed magnetic structure only in (Ti,Co)O2 anatase but not in rutile phaseMatsumoto et al. (2002), seemingly suggesting that the magnetization of the impurity atoms depends sensitively on the local environment surrounding them. More recently, Park $`et`$ $`al`$. have successfully grown ferromagnetic (Ti,Co)O2 rutile films by sputteringPark et al. (2002b). The Curie temperature was estimated to be above 400 K for Co content of 12$`\%`$. Since the rutile phase is thermodynamically more stable than the anatase phase, it may have higher potential for technical applications. The dissimilar behaviors of films grown by different methods indicate that the local environment of Co could have drastic effect on the magnetic structure of this material. Furthermore, by annealing the (Ti,Co)O2 samples prepared by the sol-gel method at different temperatures, Soo $`et`$ $`al`$. found that the local structure arround Co could remain virtually anatase-like while the nanocrystalline materials underwent anatase-to-rutile phase transitionSoo et al. (2002). This phenomenon is rather unique and certainly deserves more detailed investigation in order to gain a microscopic understanding of the magnetism displayed by this intriguing material.

We here point out that in extracting the information of the local structure of Co through the X-ray absorption fine structure method, Soo $`et`$ $`al`$. employed a refernce bulk system assuming all Co ions substitute for the Ti sites in rutile TiO2Soo et al. (2002). In a recent work, we showed by $`ab`$ $`initio`$ density functional theory study that Co could also occupy the interstitial sites in anatase TiO2. The question arise is: Is it also the case for Co-doped rutile TiO2? To answer this question, we have performed a computational investigation on the structural, electronic, and magnetic properties of Co-doped TiO2 rutile. We show that the substitutional Co ions tend to cluster, but this has no remarkable effect on the magnetization of Co. We also find that an interstitial site will become energetically favorable for a Co atom (in reference to bulk Co) when two of its neighboring substitutional sites are occupied by Co. The interstitial Co is in the low spin state (0.94$`\mu _B`$), but it exerts strong detrimental effect on the magnetization of its neighboring Co (0.76$`\mu _B`$ $``$ $`0.03\mu _B`$) and therefore reduces the magnetic moment Co in this magnetic semiconductor. In the following, we give a brief outline of our theoretical procedure, a discussion of the results based on bonding characteristics and its relevance to our understanding of the magnetism in this Co-doped oxide.

The calculations were done using a spin-polarized version of the Vienna ab initio Simulation Package (VASP)Kresse and Furthmuller (1996). Fully nonlocal Vanderbilt-type ultrasoft pseudopotentialsVanderbilt (1990) were used to represent the electron-ion interaction. Exchange correlation interactions were described by the Perdew-Wang 1991 generalized gradient approximation (GGA)Perdew et al. (1990). Wave functions and augmentation charges were expanded in terms of plane waves with an energy cutoff of 396 and 700 eV, respectively. For each system, the geometry was relaxed until the atomic forces were smaller than 0.03 eV/Å.

The 6-atom rutile unit cell of TiO2 (see Figure 1a) is characterized by the two lattice constants $`a`$ and $`c`$ and the internal parameter $`u`$. The positions of Ti are (0, 0, 0) and (1/2, 1/2, 1/2); and those of O are ($`\pm u`$, $`\pm u`$, 0) and (1/2$`\pm u`$, 1/2$`u`$, 1/2). The optimized $`a`$, $`c`$, and $`u`$ are 4.661 Å, 2.973 Å, and 0.638, in good agreement with experimentWyckoff (1960) and previous pseudopotential GGA calculationsGoniakowski et al. (1996). For Co-doped cases, we mainly dealt with a supercell containing 16 primitive unit cells ($`2a\times 2a\times 2c`$, see Figure 1b) with one or two Ti atoms (site B1, or, sites B1 and B2) replaced by Co and with/without an interstitial Co (site A). Base on our knowledge of anatase (Ti,Co)O2 that volume contraction and expansion associated with substitutional and interstitial Co have only negligible effect on the formation heat of (Ti,Co)O2Geng and Kim (2002), we fixed $`a`$ and $`c`$ at the calculated bulk TiO2 values and optimized only internal freedoms upon Co-doping.

To see whether substitutional Co (CoS) cluster or not, we first minimized the free energy of the cell Ti7CoO16 (Co content: 12.5$`\%`$). Then we doubled this cell and brought the two CoS together to form Ti14Co$`{}_{B1}{}^{}{}_{}{}^{S}`$Co$`{}_{B2}{}^{}{}_{}{}^{S}`$O32 ($`2a\times 2a\times 2c`$). We find the latter configuration is more stable by 0.31eV/Co, indicating that Co clustering does occur at this concentration. The calculated formation energy of an interstitial Co (CoI) as a neutral defect in reference to bulk Co, $`E`$(Ti16Co$`{}_{A}{}^{}{}_{}{}^{I}`$O32)-$`E`$(Ti16O32)-$`E`$(Co), is +1.75 eV, implying a sole CoI is highly unstable in rutile TiO2. This can be understood in the bond-order$``$band-strength picture. Ti:O + Co:O is much weaker than O:Ti:O, as is evident from the comparison of the formation heat of TiO2, TiO, and CoOCRC (1995). Nonetheless, with the presence of Co$`{}_{B1}{}^{}{}_{}{}^{S}`$, the formation energy of CoI, i.e., $`E`$(Ti15Co$`{}_{B1}{}^{}{}_{}{}^{S}`$Co$`{}_{A}{}^{}{}_{}{}^{I}`$O32)-$`E`$(Ti15Co$`{}_{B1}{}^{}{}_{}{}^{S}`$O32)-$`E`$(Co), drops to +0.47 eV. According to our calculations, the formation heat of CoO2 in a hypothetical rutile TiO2 structure is $`3.46`$ eV, slightly higher than two times of that of CoO ($`1.59`$ eV). This means that Co:O + Co:O is only slightly weaker than O:Co:O , thus CoI becomes much less unstable. Further calculations show that when both B1 and B2 site are occupied by Co upon clustering, site A becomes favorable for a CoI with a formation energy of $`0.72`$ eV. Once again, this can be explained by the higher formation heat of Co3O4 ($`9.57`$ eV) than two times of that of CoO2 ($`3.46`$ eV). Our calculations thus reveal a 2CoS+CoI local structure in Co-doped TiO2 rutile. It should be noted that the formation of more complicated (CoS, CoI) combinations is also possible depending on the growth condition, but the formation probability should be lower than 2CoS+CoI.

In order to examine the effect of clustering on magnetism, we performed calculations on both intra-cell FM and intra-cell AFM alignments of the Co pair in a Ti14Co$`{}_{B1}{}^{}{}_{}{}^{S}`$Co$`{}_{B2}{}^{}{}_{}{}^{S}`$O32 unit cell (For simplicity, we denote this cell as Ti14Co2O32 in the following) while keeping the inter-cell coupling ferromagnetic. The AFM phase is found to be unstable and it switches to FM during self-consistency iterations. When the two Co ions are far apart by 7.23 Åwith one at (0, 0, 0) and the other at ($`a`$, $`a`$, $`c`$), the FM phase is lower in totl energy by 0.26 eV/Co than the AFM phase. In the case of 2CoS+CoI, we checked all the three possible intra-cell spin alignments of the three Co ions. All of them converge to the same final state, with both CoI and CoS in the low spin state. The calculated total spin magnetic moment in a unit cell of Ti15Co1O32, Ti14Co2O32, and Ti14Co3O32, and that in the atomic sphere of Co (radius=1.10 Å) are listed in Table I. It is seen that in all three cases, Co ions are in the low spin state, regardless to their oxidation states. The spin moment of CoS does not change much upon clustering. CoI, on the other hand, behaves rather differently. When coupled with CoS, CoI kills almost entirely the spin moments of CoS. As a result, the average spin magnetic moment of Co is reduced to 0.29 $`\mu _B`$, less than a half of the non-CoI case.

Figure 2 displays the calculated density of states for non-doped and various cases of Co-doped rutile. Like in anatasePark et al. (2002a); Geng and Kim (2002), the Co states are mainly located in the energy-gap region and the O and Ti states are not much affected by Co doping. It is worth noting that with GGA, while Ti15Co1O32 and Ti14Co2O32 just miss being insulators, the Fermi energy of their rutile counterparts just overpass the valence bands and fall into the energy gap. The presence of CoI does not change the insulating electronic structure of this material although sole CoI ions would turn the material into a metallic system (DOS not shown), and as a consequence the magnetic interaction between Co ions at low temperature remains to be superexchange in the absence of carriers.

In Figure 3, we compare the $`3d`$ DOS of Co in rutile with with those of Co in other Co oxides with different oxidation states such as CoO2 (+4) and Co3O4 (+3 and +2). CoO2, the end member of LiCoO2, has recently been isolatedAmatucci et al. (1996). It seems to be of CdCl2 type with a monoclinic distortion, but positions of O remain to be determined. In this regard, we believe a study on a hypothetical CoO2 in TiO2 rutile structure is meaningful in illustrating the electronic structure of Co 4+. For both CoO2 and Co3O4, we have optimized their internal structural parameters. We can see that although much more localized, $`d`$ DOS of CoS in Ti15Co1O32 resembles that of Co in CoO2 in both exchange splitting and crystal field splitting, suggesting they are in +4 oxidation state. In an octahedral field, $`t_{2g}`$ is lower than $`e_g`$ and thus Co4+ will have a spin magnetic moment of 1$`\mu _B`$. As for the combination 2CoS+CoI, a close resemblance can be found between CoS (CoI) in TiO2 rutile and CoB (CoA) in Co3O4. It thus provides a strong evidence that CoS and CoI in such a combination are in +3 and +2 oxidation states, respectively. The six $`d`$ electrons of CoS just fill $`t_{2g}`$ states and thus yield a vanishing spin magnetic moment. Unlike the interstial Co in anatase and CoA in Co3O4, CoI in TiO2 rutile is in an octahedral field, therefore the lone electron in $`e_g`$ state gives a spin moment of 1$`\mu _B`$.

In summary, our $`ab`$ $`initio`$ density functional theory study reveals that in Co-doped rutile TiO2, Co can occupy interstitial sites, upon clustering of sbstitutional Co. As a result, a reference bulk system used to extract information of the local structure of Co in a polycrystalline material should be one that containing both substitutional and interstitial Co. It is found that the interstitial Co destroys the spin-polarization of the neighboring substitutional impurities by lowering its charge state and therefore plays an important role in the magnetic structure of such material.

###### Acknowledgements.

Work in Korea was supported by Korea Institute of Science and Technology Evaluation and Planning (Creative Research Initiative) and in Germany by Max-Planck-Society Fellowship. W.T.G. thanks Dr. Qiang Sun for helpful discussions.


text:
# Influence of progressive central hypovolemia on multifractal dimension of cardiac interbeat intervals.

## I Introduction

In the past two decades we have witnessed an explosion in the biophysics and physiological literature with regard to the identification of phenomena having long-term memory and probability densities that extend far beyond the typical tail region of Gaussian distributions. One way these processes have been classified is as 1/f-phenomena, since their time series have spectra that are inverse power law in frequency or their probabilities have an inverse power-law distribution. In either case the underlying structure is fractal, either in space, time or both bass94 . Physiological signals, such as the heartbeat interval time series studied herein, are typically generated by complex self-regulatory systems that process inputs with a broad range of characteristics. Such physiological time series fluctuate in an irregular and complex manner such as shown in Fig. 1 and the statistics of these fluctuations often exhibit self-affine or fractal properties bass94 . Such processes do not possess dominant time scales, but they can frequently be characterized by fractal dimensions feders . Ivanov et al ivanov99  established that healthy human heartbeat intervals, rather than being fractal, exhibit multifractal properties and they uncovered the loss of multifractality for a life-threatening condition known as congestive heart failure. Latka et al. latka02  similarly determined that cerebral blood flow in healthy humans is also multifractal.

Recently, Scafetta et al. scafetta1  studied how the fractal and multifractal nature of the human stride interval is modified by changing the pace velocity in both free and metronomically constrained walking conditions and introduced a neuronal super central pattern generator (SCPG) to model the human stride interval time series scafetta2 . In particular, it was determined that an increase of neural stress, induced by a faster or slower pace than normal, causes a slight increase of both fractal and multifractal properties of the free pace stride interval time series scafetta1 ; scafetta2 . In the present analysis we used lower body negative pressure (LBNP) to induce progressive central hypovolemia and determined how such stimulation of autonomic neural control mechanisms modified the statistical nature of heart-beat time series. A LBNP is applied by sealing supine subjects at the level of the iliac crest in an air-tight chamber connected to a vacuum pressure control system. The magnitude of LBNP is carefully regulated to induce progressive central hypovolemia and autonomic cardiovascular compensatory mechanisms are evaluated. LBNP decreases central venous pressure, stroke volume and cardiac output by inducing a footward fluid shift, and increases heart rate and total peripheral resistance through parasympathetic neural withdrawal and sympathetic neural activation Covertino . The magnitude of changes in autonomic neural activity can be estimated non-invasively in the frequency domain, but common analyses incorporating either fast-Fourier based power spectra Akselrod  or autoregressive modeling Malliani  do not account for inherent multifractality or non-stationarity underlying complex cardiovascular regulatory mechanisms. For this reason, in the present study we applied two wavelet-based techniques to study the fractal and multifractal dependence of the cardiovascular responses to progressive central hypovolemia to determine how cardiovascular stress affects fluctuations of the data on a hierarchy of time scales. For example, the analysis allows us to isolate the influence of the breathing cycle component on heart rate. We determine that in this situation, similar to the situation in gait dynamics scafetta1 ; scafetta2 , an induced stress, in this case through central hypovolemia, causes a slight increase in both the fractal and multifractal properties of the heartbeat interval time series. We suggest that such changes in fractal and multifractal properties characterize neural network responses of biological systems to physiological stress.

In Sec. II we present the heart interbeat interval data analyzed herein; Sec. III is devoted to a review of the techniques of analysis; Sec. IV is devoted to a detailed analysis of the heartbeat interval time series of 12 normal subjects and in Sec. V we draw some conclusions.

## II Heart interbeat data

In the main experiment we studied twelve healthy men between the ages of 26 and 40 after they had been informed of the nature of the experiment and had signed a consent form, approved by the institutional human use committee of the U.S. Army Institute of Surgical Research, Ft Sam Houston, TX. Subjects were positioned supine in an LBNP chamber, sealed at the waist to the level of the iliac crest, and instrumented with a standard 4-lead electrocardiogram. Because respiratory frequency and depth may seriously confound interpretation of autonomic rhythms Brown , subjects breathed in synchrony with a metronome set at a pace of 15 breaths per minute (0.25 Hz). Electrocardiogram data were sampled at 500 Hz during 3-minute periods of controlled frequency breathing at negative pressure of 0, -15, -30, -45, -60, and again at 0 mmHg. Data were recorded directly to computer (Windaq, Dataq Instruments, Akron, OH) and then imported into a customized software program for analysis \[Cardio-Pulmonary Research Software for Windows, (WinCPRS), Absolute Aliens Ay, Turku, Finland\]. R-waves were automatically detected from continuous waveforms and marked at their occurrences in time to obtain heartbeat periods (tachogram) as shown in Fig. 1.

In a second experiment, conducted at the Institute of Clinical Physiology, CNR-Pisa, we compared the fluctuations of a tachogram with the contemporary respiratory signal regarding one person. This experiment was originally done with thirteen healthy volunteers. These volunteers underwent an experimental session in which an electrocardiographic (ECG) derivation and the respiratory signal were simultaneously and continuously recorded at 1000 Hz sampling frequency. The ECG was recorded using standard electrodes, while the respiratory signal was detected through a polymeric piezoelectric dc-coupled transducer inserted into a belt wrapped around the chest. During the experimental session the subjects were comfortably sitting and breathing freely.

## III Fractal and multiresolution analysis based on wavelets

The analysis we apply to heartbeat time series emphasizes the scaling nature of the phenomenon. In the next two subsections we briefly review the two techniques we use in our analysis. The first method determines the distribution of Hölder exponents of a time series, by means of the continuous wavelet transform (CWT) struzik . This technique is used to study the fractal and multifractal nature of a time series and has already been used to analyze the human stride interval time series scafetta1 . The second method uses the wavelet multiresolution analysis (WMA) of a time series, by implementing the maximum overlap discrete wavelet transform (MODWT) percival . The latter technique provides a detailed study of the fluctuations of a time series on a scale-by-scale basis.

### III.1 Approximate estimation of local Hölder exponents and their probability distribution.

Mandelbrot 2Mandelbrot  showed that many natural phenomena are described by self-affine, correlated, scaling time series and as we mentioned above, heartbeat data fall into this category. The scaling properties of the fractal noise studied by Mandelbort, fractional Gaussian noises (fGn) 2Mandelbrot , are characterized by an exponent that he called $`H`$ in honor of Hurst. Consequently, if $`X(t)`$ is a fractal process with Hurst exponent $`H`$ and $`c`$ is a constant, then $`X\left(t\right)=X(ct)/c^H`$ is another fractal process with the same statistics. The phenomenon of $`fGn`$ has a spectrum with the inverse power-law form

$$P(f)f^\beta =f^{12H}f^{12h_0},$$
(1)
where $`f`$ is the frequency, $`H`$ is the Hurst exponent and $`h_0`$ is the average of the Hölder exponent distribution among the singularities of the time series. We can also express a fractal noise $`\{\xi _i\}`$ through the equivalent representation of the autocorrelation function politi

$$C(r)=\frac{\xi _i\xi _{i+r}}{\xi _i^2}r^{2H2}r^{2h_0}.$$
(2)
The self-affine property expressed by (1) and the relation between $`\beta `$, $`H`$ and $`h_0`$ are theoretically valid only for a infinitely long monofractal time series. Figure 2 shows a computer generated realization of $`fGn`$ with Hurst exponent $`H=1`$ or Hölder exponent $`h_00`$, also known as $`1/f`$ noise or pink noise. This type of noise is important because it represents a threshold between the persistent-stationary noise ($`0.5<h_0<0`$) and the non-stationary noise ($`h_0>0`$). Moreover, the $`1/f`$ noise is characterized by an autocorrelation function $`C(r)`$ that does not depend on the distance $`r`$, Eq. (2). The random noise or white noise is characterized by $`h_00.5`$ whereas the random walk or Brownian motion is characterized by $`h_00.5`$. In general, the higher the Hölder exponent, the smoother the time series.

As seen in Figure 2, a fractal noise is characterized by trends and discontinuities that give a particular geometric shape to the signal. The rapid changes in the time series are called singularities of the signal and their strength is measured by a Hölder exponent struzik . Given a function $`f(x)`$ with a singularity at $`x_0`$, the Hölder exponent $`h(x_0)`$ at such a point is defined as the supremum of all exponents $`h`$ that fulfills the condition:

$$|f(x)P_n(xx_0)|C|xx_0|^h,$$
(3)
where $`P_n(xx_0)`$ is a polynomial of degree $`n<h`$.

The Hölder exponent of a singularity can be evaluated by using the wavelet transform daubechies ; Mallat ; percival . Wavelet transforms make use of scaling functions that have the property of being localized in both time and frequency. A scaling coefficient $`s`$ characterizes and measures the width of a wavelet. Given a signal $`f(x)`$, the continuous wavelet transform (CWT) of $`f(x)`$ is defined by

$$W_{s,x_0}(f)=\underset{\mathrm{}}{\overset{\mathrm{}}{}}\frac{1}{s}\psi \left(\frac{xx_0}{s}\right)f(x)𝑑x,$$
(4)
where the kernel $`\psi (u)`$ is the wavelet filter centered at the origin, $`u=0`$, with unit width, $`s=1`$. A widely adopted choice of the kernel $`\psi (u)`$ is the so-called Mexican Hat, that is, the second derivative of a Gaussian.

The wavelet transform can be used to determine the Hölder exponent of a singularity because the wavelet kernel $`\psi (u)`$ can be chosen in such a way as to be orthogonal to polynomials up to degree $`n`$, that is, such that the following properties are fulfilled

$$\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\frac{1}{s}\psi \left(\frac{xx_0}{s}\right)x^m𝑑x=0m,0mn.$$
(5)
In fact, if (5) holds true, it is easy to prove that if the function $`f(x)`$ fulfils condition (3), its wavelet transform at $`x=x_0`$ is given by

$$W_{s,x_0}(f)=C|s|^{h(x_0)}\underset{\mathrm{}}{\overset{+\mathrm{}}{}}|u|^{h(x_0)}\psi (u)𝑑u|s|^{h(x_0)},$$
(6)
where $`u=(xx_0)/s`$. Therefore, at least theoretically, the Hölder exponent of a singularity that is localized in $`x_0`$, can be evaluated as the scaling exponent of the wavelet transform coefficient, $`W_{s,x_0}(f)`$, for $`s0`$.

Even if Eq. (6) can be evaluated for any position $`x_0`$, we are interested only in the cusp singularities of the time series. Mallat et al. Mallat ; mallat2 ; mallat3  show that the Hölder exponent of these singularities can be evaluated by studying the scaling exponent $`h(x_0)`$ along the so-called maxima line that converges towards the singularity. The maxima lines are defined by the extremes of the wavelet transform coefficients (4) at each wavelet scale $`s`$. Arneodo et. al. arneodo1 ; arneodo2  proved that WTMM can be used to define a multifractal-like formalism that gives the stochastic properties of the singularities of a fractal or multifractal noise. This methodology has been recently used to determine the multifractal nature of many signals, for example, that for human heartbeats ivanov99 .

However, the above method presents some problems of stability when applied to observational data. In fact, the spectrum can be corrupted by the divergences of negative moments Mallat ; arneodo2 ; struzik2  or by the outliers, that is, the end points of the sample singularities struzik . Different methods have been suggested to remove the divergences due to the negative moments of the multifractal partition function, for example, by chaining the wavelet maxima across scales Mallat  or, more efficiently, by bounding the Hölder exponent of the maxima line by using the slope wavelet struzik2 . The corruption of the singularity spectrum due to the outliers is more difficult to deal with.

Struzik struzik  suggested an alternative method that has the ability to determine an approximate value of local singularity strength. The spectrum may then be evaluated from these approximate values. The idea is to estimate the mean Hölder exponent $`\overline{h}`$ as a linear fit of the following equation

$$\mathrm{log}[M(s)]=\overline{h}\mathrm{log}(s)+C,$$
(7)
where the function $`M(s)`$ is obtained via the partition function, and the mean Hölder exponent $`\overline{h}`$ is the local version of the Hurst exponent $`H`$ 2Mandelbrot . More precisely, for a monofractal noise with Hurst exponent $`H`$, we have $`\overline{h}=H1`$ because the Hurst exponent is evaluated by integrating the noise 2Mandelbrot ; feders . Here again the equality is only rigorously true for an infinitely long monofractal noise data set. The approximate local Hölder exponent $`\widehat{h}(x_0,s)`$ at the singularity $`x_0`$ can now be evaluated as the slope

$$\widehat{h}(x_0,s)=\frac{\mathrm{log}(|W_{s,x_0}(f)|)(\overline{h}\mathrm{log}(s)+C)}{\mathrm{log}(s)\mathrm{log}(s_N)},$$
(8)
where $`s_N`$ is the length of the entire wavelet maxima line tree, that is, the maximum available scale that coincides with the sample length $`s_N=N`$, and $`x_0`$ belongs to the set $`\mathrm{\Omega }(s)`$ of all wavelet maxima at the scale $`s`$ that assume the value $`W_{s,x_0}(f)`$.

While further details may be found in the papers by Struzik struzik  and in Ref. scafetta1 , here we focus on the interpretation of the final output, that is, the distribution of estimated local Hölder exponents of a time series as given in scafetta1 . Fig. 3 shows the histogram of Hölder exponents obtained using a computer-generated data set with a Hurst exponent $`H=1`$ shown in Fig. 2. The histogram is fitted with a normalized Gaussian function of the type

$$g(h)=\frac{1}{\sqrt{2\pi }\sigma }\mathrm{exp}\left[\frac{(hh_0)^2}{2\sigma ^2}\right]$$
(9)
that captures the two main characteristics of the time series: (a) the average of Hölder exponent $`h_0`$ that gives an estimate of the fractal nature of a time series; (b) the width of the distribution given by the standard deviation $`\sigma `$ that measures the variability of the local Hölder exponents, that is, the multifractality of a time series. Finally, as Scafetta et al. scafetta1  explained, a monofractal time series presents a particular width of the Hölder exponent distribution that depends on the length of the time series, as Eq. (8) suggests. Therefore, a time series will present a multifractal nature only if the width of its Hölder exponent distribution is larger than that of a monofractal time series of the same length. In general, given two time series of the same length, if $`\sigma _1>\sigma _2`$, we may say that the time series ‘$`1`$’ has a stronger multifractal character than time series ‘$`2`$’.

### III.2 Wavelet multiresolution analysis.

The second wavelet-based technique that we implement allows a detailed study of the fluctuations of a time series on a scale-by-scale temporal basis. The MODWT percival  is almost independent of the particular family of wavelets used in the analysis and is the basic tool needed for studying the multiresolution analysis of time series of $`N`$ data points via wavelets. In the book of Percival and Walden, percival , the reader can find all the mathematical details.

In this paper we make use of the Daubechies least asymmetric scaling wavelet filter (LA8) that looks like the Mexican Hat, but is also weakly asymmetric; a fact that makes LA8 filter more malleable than the Mexican Hat. The WMA via MODWT establishes that given an integer $`J`$ such that $`2^J<N`$, where $`N`$ is the number of data points, the original time series represented by the vector $`𝐗`$ can be decomposed on a hierarchy of time scales represented by a smooth part plus details as follows:

$$𝐗=S_J+\underset{j=1}{\overset{J}{}}D_j,$$
(10)
with the quantity $`S_j`$ generated by the recursion relation

$$S_{j1}=S_j+D_j.$$
(11)
The detail $`D_j`$ of Eq. (10) represents changes on a scale of $`\tau =2^j`$, while the smooth $`S_J`$ represents the smooth wavelet averages on a scale of $`\tau _J=2^J`$. In the same way we refer to as residuals the curves

$$R_J=\underset{j=1}{\overset{J}{}}D_j,$$
(12)
such that $`𝐗=S_J+R_J`$.

Fig. 4 shows that WMA of the heartbeat data shown in Fig. 1 for J=4. The detail curves show the main characteristic fluctuation that characterize each time scale that may be expressed in the number of beats, as done in Fig. 4, or in physical time units. We shall use these detail curves to establish the dependence on the pressure of the period measured in both number of beats and in physical time units of the fluctuations of the heartbeat at different wavelet scale.

## IV Analysis of the heartbeat time series

Fig. 1 shows the heartbeat time series at the four conditions of LBNP for a typical individual. In this work, we analyze the heartbeat interval time series for 12 individuals. In the time series data shown in the figure the increase in the average heart rate with increasing negative pressure is clear, as is the return to normal with the return of the negative pressure to zero. To emphasize this effect, in Fig. 5 we plot the average heartbeat interval against the pressure for 12 individuals. The figure shows that the average heartbeat interval manifests a steady decrease from $`\overline{\mathrm{\Delta }t}980\pm 50`$ msec at $`P=0`$ mmHg to $`\overline{\mathrm{\Delta }t}710\pm 50`$ msec at $`P=60`$ mmHg. At the recovery, the average heartbeat interval is slightly larger than the initial one; we measure $`\overline{\mathrm{\Delta }t}1050\pm 50`$ msec. The error bars measure the average standard deviation of the variability of the heartbeat interval for each individual and do not seem to significantly change with the pressure.

A careful visual analysis of Fig. 1 suggests that by increasing the negative pressure, the heartbeat interval time series becomes smoother. As explained in the previous section, a decrease of randomness or an increase of the smoothness of a time series may be easily detected by the Hölder exponents of the time series that measure the strength of the singularities. The Hölder exponents increase with the smoothness of a time series. Fig. 6 shows an estimation of the Hölder exponents of the time series of the heartbeat interval time series shown in Fig. 1. We use the algorithm discussed in the previous section to obtain these exponents. Fig. 6 shows that the Hölder exponents significantly increase with the pressure, indicating that the fractal dimension of the heartbeat interval monotonically decreases with negative pressure. At recovery, the fractal properties return to their original values.

Fig. 7 depicts the average Hölder exponent, indicated by the mean value $`h_0`$, and by the standard deviation $`\sigma `$ for the 12 individuals against the pressure. As explained in the previous section the standard deviation measures the variability of the Hölder exponent values that can be related to the multifractality of a time series. The data indicates that the average mean value has a steady increase from $`\overline{h_0}0.177\pm 0.072`$ at $`P=0`$ mmHg to $`\overline{h_0}0.332\pm 0.104`$ at $`P=60`$ mmHg. The standard deviation is indicated by the error bars. At the recovery when the negative pressure returns to zero we get $`\overline{h_0}0.157\pm 0.062`$. We also observe that by increasing the negative pressure the average standard deviation tends to increase. The values we measure are $`\overline{\sigma _0}=0.072`$, $`\overline{\sigma _{15}}=0.062`$, $`\overline{\sigma _{30}}=0.066`$, $`\overline{\sigma _{45}}=0.077`$, $`\overline{\sigma _{60}}=0.104`$ and $`\overline{\sigma _0}=0.062`$ at the recovery. As anticipated, an increase of the standard deviation of the Hölder exponents for equal length time series may be interpreted as an increase of multifractal properties of the time series. Even if the time series taken for each value of the pressure are relatively short (approx. 200 data points for $`P=0`$ mmHg and approx. 300 data points for $`P=60`$ mmHg) and, therefore, the statistics may be not optimum. However, the almost monotonic increase of the standard deviation with the pressure followed by an abrupt decrease at the recovery may reasonably indicate a reliable increase of the multifractality of the heartbeat interval time series with increasing negative pressure.

To emphasize the properties of heartbeat data at various temporal resolutions we need to make use of WMA of the data. In Fig. 4 we show the WMA of the heartbeat time interval depicted in Fig. 1. These data are analyzed using the MODWT with the LA8 wavelet filter for four levels of resolution, that is, until the wavelet scale $`J=4`$. As explained in the previous section the smooth curve S4 captures the smooth trend of the fluctuation at the scale of $`\tau =2^4=16`$ consecutive heartbeats. The details curves D1-D4 capture the local fluctuations of the signal at each of the four time scale $`\tau =2^i`$. This means, for example, that D1 captures the fluctuations of the original signal with a period that is approximately between 2 and 4 beats. We recall that the shortest possible period that can be detected in a time series is 2 time units. The detail D2 captures the fluctuations with a period in the interval between 4 and 8 beats, and similarly for the following details. According to Eq. (10) the sum of the four details D1-D4 and the smooth curve S4 gives the original datasets shown in Fig. 1.

Fig. 4 stresses the great utility of analyzing times series of complex systems through the WMA via the MODWT. In fact, the wavelet sensibility to the local changes of the signal allows us to easily extract information that is impossible to obtain by using, for example, the Fourier transform since the latter averages the changes of the entire signal at each Fourier transform frequency. For example, Fig. 4 clearly shows that the amplitude of the local decomposition is not constant but fluctuates in time. In particular we notice that the amplitude of the fluctuations captured by the details D3 and D4 does not seem to change systematically with the pressure. Therefore, we can conclude that the change of pressure do not influence significantly the slow fluctuations with a period larger than 8 heartbeats.

On the contrary, the details D1 and D2 clearly show a significant reduction of the amplitude of the fast fluctuation with increasing negative pressure. At the recovery, the amplitude of the fluctuations returns to a value similar to the initial one for $`P=0`$ mmHg. The reduction of the amplitude of the fast fluctuations manifests itself with the increase of the smoothness of the heartbeat interval time series that we have already noticed to occur by increasing the negative pressure, as Fig 1 shows. We have also noticed that an increase of the smoothness yields an average increase of the local Hölder exponents, that are associated with the fractal dimension of a time series, as explained earlier and shown in Figs. 6 and 7. Fig. 4 shows the importance of using WMA because we are able to isolate the component of the signal of the heartbeat interval time series that presents dependence on the pressure and identifies it with the changes captured by the wavelet details D1 and D2.

Because each detail curve captures the variation of the signal inside a frequency or period range, it is possible to also study how the average period of the fluctuations at a particular wavelet scale are affected by the pressure. Fig. 8 shows the average mean period measured in number of beats of the fluctuation associated to each detail curve obtained by the WMA for 12 subjects. We stress that these calculations are made without any kind of interpolation between the data. So the natural unit of the time sequence is the number of heartbeats. To get the results in physical time units we have to convert the results obtained with the heartbeat number in physical time units by using the original R-R sequence or we may attempt a direct calculation of WMA after an interpolation procedure that generates an even time series in a physical units. However, interpolation alters data, so we prefer the first procedure. Fig. 9 shows the average mean period measured in physical time units of the same fluctuations. We decided to plot the results concerning D1, R2, D3 and D4 curves. We stress that we are considering the residual R2 that is given by $`R_2=D_1+D_2`$ according to (11) rather than the D2 curve. We prefer to use the residual R2 instead of the detail D2 because the range from 2 to 8 heartbeats is supposed to capture the high frequency component of the R-R signal that is coupled with the respiratory signal. The residual R2 fully contains such information that, instead, is split between D1 and D2 curves with D1 being a kind of harmonic of D2. For this experiment we recall that the subjects breathed in synchrony with a metronome sets a pace of 15 breaths per minute (0.25 Hz). So we expect to find such a frequency in the R2 curves. The detail D3 captures the range from 8 to 16 heartbeats and corresponds to the low frequency component of the R-R signal. Higher levels of details, that are included in the S3 smooth curve, capture the very-low component of the interbeat signal.

The results plotted in Figs. 8 and 9 concerning the detail curves D3 and D4 show a more erratic behavior with large errorbars and their interpretation remains unclear. Fig. 8 shows an increase of the heartbeat number for cycles captured by D3 curve with increasing negative pressure. Fig. 9 shows that in physical time unit the period of these D3 cycles remains almost constant at almost 10 seconds, that is, at a frequency of 0.1 Hz, while the period of the cycles associated with D4 detail curve tends to decrease with increasing negative pressure.

The results concerning the fast fluctuations, curves D1 and R2, show more regular patterns. In particular we notice that by increasing the negative pressure, Fig. 8 shows that the average mean period of the fluctuations measured in heartbeat number is almost constant for details D1, D3 and D4. Instead, the results concerning the R2 residual shows a sensitive steady increase of the average number of heartbeats from almost 4.1 heartbeats per cycle to 4.9 heartbeats per cycle with increasing negative pressure from $`P=0`$ mmHg to $`P=60`$ mmHg. At recovery the average number of heartbeats returns to 4.1 heartbeat per cycle. By looking at the same results expressed in physical time units, see Fig. 9, we observe that the average period of the fluctuations captured by the R2 residual curve is almost constant with the pressure. We measure 4 seconds that correspond to a frequency of 0.25 Hz. Instead, the average period of the fluctuations captured by the D1 detail curve shows a steady decrease from 3.2 to almost 2.1 seconds with increasing negative pressure. The different behaviors characterizing the wavelet detail D1 and R2 depicted in Figs. 8 and 9, suggest that the two curves can be partially decoupled because they respond differently to the increases in negative pressure.

Fig. 10 shows a comparison between the R2 residual curve of the heartbeat interval and respiratory signals of a person for 60 seconds under normal conditions. These data were obtained in a second experiment. Figure 10 shows a very strong correspondence between the fluctuations captured by the R2 residual curve and the respiratory cycle. This comparison suggests that the R2 residual curve has the physiological meaning of describing the fluctuations of respiratory sinus arrhythmia. In fact, as Fig. 9 shows, the fluctuations captured by the R2 residual have the frequency of 0.25 Hz of the respiratory signal that is kept constant during the hypovolemia experiment.

## V Summary and conclusions

The analysis of the data suggests that progressive central hypovolemia induced with LBNP:

1) Causes decreases of R-R intervals and, therefore, increases of heartbeat rate.

2) Decreases heart rate variability and therefore smoothes out the heartbeat time series.

3) Reveals that two measures of the decrease in the heartbeat rate variability with increasing negative pressure are an increase in the Hölder exponent and a decrease in the fractal dimension of the time series.

4) Reveals that accompanying the increase in the Hölder exponent, with increasing negative pressure, is an increase in the multifractality of the time series as measured by the width of the Hölder exponent distribution. This non-stationary aspect of the data was also observed in the study of human gait data, when the walkers were put under physical stress by asking them to walk faster or slower than normal scafetta1 ; scafetta2 .

5) Reveals that independent measures of the heartbeat rate response to the increases in negative pressure are the wavelet transform details D1 and D2. Both these non-periodic functions display a monotonic decrease in amplitude with increasing negative pressure.

6) Suggests that the amplitude of the higher-order details, D3 and D4, are independent of the changes in the negative pressure. However, the period of the cycles captured by D4 tends to decrease with increasing negative pressure. This suggests that part of the spectral energy of the very-low frequency ($`f<0.04`$ Hz) component of the RR signal moves toward the low frequency ($`0.04<f<0.15`$ Hz) component of the same signal.

7) Suggests that the wavelet transform residual curve R2 is driven by the respiratory cycle. The residual R2 captures the high frequency component of the RR signal ( $`0.15<f<0.4`$ Hz). This is an indication that the wavelet transform is able to sort through the fluctuations of the cardiac time series and highlight particular physiological phenomena.

These separate observations made from the analysis of the heartbeat time series stressed by central hypovolemia are interesting in themselves, but they suggest a much more significant finding. They suggest that the influence of stress on a physiological system is manifest in the non-stationary character of the time series and that the degree of non-stationarity can be measured through the level of fractality and multifractality of the time series regarding the stressed system.

Acknowledgment: This work was supported in part by a Faculty Research Fellowship through the US Army Institute of Surgical Research under the auspices of the US Army Research Office Scientific Services Program administered by Battelle (DAAH04-96-8806). We acknowledge Dr. Victor Convertino, Institute of Surgical Research, Ft. Sam Houston TX, as having contributed the LBNP data. We thank the research group of prof. B Ghelarducci, Dipartimento di Fisiologia e Biochimica, Universitá di Pisa, and Dr. M. Varanini, Institute of Clinical Physiology, CNR-Pisa, for the respiratory and contemporary heartbeat signal used in the second experiment. N.S. thanks the Army Research Office for support under grant DAAG5598D0002.