extfieldsofstudy:
Biology
Computer Science
Medicine

text:
A configuration space of homologous proteins conserving mutual information and allowing a phylogeny inference based on pair-wise Z-score probabilities

Background Popular methods to reconstruct molecular phylogenies are based on multiple sequence alignments, in which addition or removal of data may change the resulting tree topology. We have sought a representation of homologous proteins that would conserve the information of pair-wise sequence alignments, respect probabilistic properties of Z-scores (Monte Carlo methods applied to pair-wise comparisons) and be the basis for a novel method of consistent and stable phylogenetic reconstruction. Results We have built up a spatial representation of protein sequences using concepts from particle physics (configuration space) and respecting a frame of constraints deduced from pair-wise alignment score properties in information theory. The obtained configuration space of homologous proteins (CSHP) allows the representation of real and shuffled sequences, and thereupon an expression of the TULIP theorem for Z-score probabilities. Based on the CSHP, we propose a phylogeny reconstruction using Z-scores. Deduced trees, called TULIP trees, are consistent with multiple-alignment based trees. Furthermore, the TULIP tree reconstruction method provides a solution for some previously reported incongruent results, such as the apicomplexan enolase phylogeny. Conclusion The CSHP is a unified model that conserves mutual information between proteins in the way physical models conserve energy. Applications include the reconstruction of evolutionary consistent and robust trees, the topology of which is based on a spatial representation that is not reordered after addition or removal of sequences. The CSHP and its assigned phylogenetic topology, provide a powerful and easily updated representation for massive pair-wise genome comparisons based on Z-score computations.

Background
Past events that gave birth to biological entities can be ten-tatively reconstructed based on collections of descriptors traced in ancient or present-day creatures. Using genomic sequences, an estimate of the relative time separating branching events, previously supported by geological records, could be formalized using mathematical models. The use of proteins for evolutionary reconstructions was vastly explored as soon as the first amino acid sequences were made available [1][2][3][4][5][6][7][8][9]. The rich biological information contained in protein sequences stems from their being, on the one hand, translation of genes that reflect the history of genetic events to which the species has been subjected, and on the other hand, effectors of the functions constituting a living creature [10] Since protein sequences are encoded in a 20-amino acid alphabet, they are also considered to embody more information-per-site than DNA or RNA [11]; they also exhibit smaller compositional trends [12,13]. When compared, sequences that share substantial features are considered as possible homologues [14], based on the fundamental postulate that can be simply stated as "the closer in the evolution, the more alike and conversely, the more alike, probably the closer in the evolution".
As summarized by Otu and Sayood [15], the techniques of molecular phylogenetic analyses can be divided into two groups. In the first case, a matrix representing the distance between each pair of sequences is calculated and then transformed into a tree. In the second case, a tree is found that can best explain the observed sequences under evolutionary assumptions, after evaluation of the fitness of different topologies. Some of the approaches in the first category utilize distance measures [16][17][18][19] with different models of nucleotide substitution or amino acid replacement. The second category can further be divided into two groups based on the optimality criterion used in tree evaluation: parsimony [20,21] and maximum likelihood methods [22,23]. For a detailed comparison of these methods see [24].
In phylogeny inference based on distance methods, features separating related proteins are used to estimate an observed distance, also called the p-distance, the simplest measure of which is just the number of different sites between proteins. Divergence time (t), also called genetic distance or evolutionary time, is calculated from the p-distance, depending on assumptions derived from evolutionary models [11,24]. For example, the assumption that mutational events happen with equal probability at each site of any sequence leads to the molecular clock model [2]. Although widely used, it is well-known to be unrealistic and numerous corrections have been proposed to refine it [19,25,26]. By definition, the distance matrix is given as T = (t ab ) where a and b represent the homologous sequences from the analyzed dataset. Tree reconstruction algorithms are then applied to these matrices [11,24]. Eventually, phylogenetic trees corresponding to the classified sequences are statistically evaluated with bootstrap methods and, when available, calibrated using dated fossils [25,26].
Doolittle et al. [27] have proposed methods for converting amino acid alignment scores into measures of evolutionary time. Similarity between amino acids [28][29][30] provides a way to weight and score alignments [31]. In practice the optimal alignment of two sequences (a and b) is determined from the optimal score s(a,b) [25,27], computed with a dynamic programming procedure [32,33]. In aligned sequences, conservation is measured at identical sites, whereas variation is scaled at substituted sites. To estimate the variation/conservation balance, the p-distance can be given as a function of f id , the fraction of identical residues: p-distance = 1 -f id . To take into account that multiple mutations can happen at the same site, an expression of f id was proposed by Doolittle et al. [27] using s(a*,b*), the score obtained from randomized a and b sequences [34] and s id , the average score of the sequences compared with themselves [19,25,27]: To connect pair-wise alignments and phylogeny, divergence time has been approximated: introducing a Poisson correction [2] as a reasonable stochastic law relating amino acid changes and elapsed time. As mentioned earlier, adjustments and corrections of equation (2) were proposed to fit more realistically the complexity of evolution [11,25,35]. This attempt of unification helped reconstructing phylogeny of major lineages [27]. However, detailed phylogenic trees obtained from evolutionary close sequences are not satisfactory. In practice, phylogenies are reconstructed based on multiple alignments. Multiple alignment based (MAB) trees are recalculated when incremented with additional sequences; although MAB methods are usually considered accurate, numerous cases of inconsistencies (incongruence) between observed data and deduced MAB trees are recorded (see [15,36]).
Here, we re-examine the estimate of the p-distance between two homologous sequences, based on f id , as a source for geometric positioning, divergence time calculations and evolutionary reconstruction. We based our model on mathematical properties that alignment scores should respect; i) information theory [37,38] applied to sequence similarity, ii) algorithmic theory applied to alignment optimization [28] and iii) alignment probability, particularly in conformity with the TULIP theorem [39]. We used these properties as a framework of constraints to build a geometric representation of a space of probably homologous proteins and define a theoretically explicit measure of protein proximity. This unified model conserves information in the way physical models conserve matter or energy. The obtained representation of protein sequences is unaltered by adding or removing sequences. Applications include therefore the reconstruction of evolutionary consistent and robust trees, the topology of which is based on a spatial representation that is not reordered after addition or removal of sequences.

Pair-wise sequence alignment scores in information theory
Criteria to measure the variation/conservation balance between proteins should embody as much as possible the structural and functional potentiality within sequences of amino acids. In the absence of explicit physical criteria, amino acid similarity was solved empirically by measuring amino acid substitution frequencies in alignments of homologous sequences [30,40]. Given two amino acids i and j, the similarity function s(i,j) was set as: where ϖ ij is the observed frequency of substitution of i by j or j by i, and π i and π j are the frequencies of i and j in the two aligned sequences. The ϖ ij frequency is the estimate of the probability of substitution of i by j in real alignments; whereas π i π j is the estimate of the probability of substitution under the independency hypothesis. The similarity function gives a 20 × 20 similarity matrix usable to score protein sequence alignments, that can be interpreted in the information theory [37,38] according to the following proposition.

Proposition 1
Amino acid substitution matrix values are estimates of the mutual information between amino acids in the sense of Hartley [37,38]. Consequently, the optimal alignment score computed between two biological sequences is an estimate of the optimal mutual information between these sequences.

Proof
Given a probability law P that characterizes a random variable, the Hartley self-information h is defined as the amount of information one gains when an event i occurred, or equivalently the amount of uncertainty one loses after learning that i happened: The less likely an event i, the more we learn about the system when i happens. The mutual information I between two events, is the reduction of the uncertainty of one event i due to the knowledge of the other j: Mutual information is symmetrical, i.e. I j→i = I i→j , and in the following will be expressed by I(i;j). The self and mutual information of two events i and j are related: If the occurrence of one of the two events makes the second impossible, then the mutual information is equal to -∞. If the two events are fully independent, mutual information is null. The empirical measure of the similarity between two amino acids described in equation (3) can therefore be expressed in probabilistic terms: where P ϖ is the joint probability to have i and j aligned in a given alignment and P π the measure of probability that amino acids occur in a given sequence. From equations (4) and (6), equation (7) becomes: that is s(i, j) = I(i; j) (9) As a consequence, the similarity function (or score) is the mutual information between two amino acids. Additionally, the score between sequences (the sum of elementary scores between amino acids, [32, 33,41,42]) is, according to the hypothesis of independence of amino acid positions, the estimated mutual information between the two given biological sequences.
Once two sequences are aligned, we pose the question whether the alignment score is sufficient to assess that the proteins are conceivably alike and thus evolutionarily related? The theorem of the upper limit of a sequence alignment score probability (TULIP theorem, [39]), sets the upper bound of an alignment score probability, under a hypothesis less restrictive than the Karlin-Altschul model [43]. Given two real sequences a and b (a = a 1 a 2 ...a m and b = b 1 b 2 ...b n ), where s = s(a,b) the maximal score of a pair-wise alignment obtained with any alignment method, b* the variable corresponding to the shuffled sequences from b, and given P{S(a,b*)≥s} the probability that an alignment by chance between a and b* has a higher score than s, then whatever the distribution of the random variable S(a,b*) the TULIP theorem states: with k > 1, µ the mean of and σ its standard deviation. The unique restriction on S(a,b*) is that it has a finite mean and a finite variance. A first corollary of the TULIP theorem is that the Z-score is a statistical test for the probability of a sequence alignment score. We additionally state the following new corollary.

TULIP corollary 2
Given the TULIP theorem conditions, let be the Z-score [44]. Then, z(a,b*) is the greatest possible value for k (k∈]1,+∞[), which holds relation (10) true. In consequence, with k = z(a,b*), then The best upper bound value for P{S(a,b*)≥s} is termed . From the TULIP theorem and corollaries, the comparison of a protein to a given reference a, weighed by an alignment score, is characterized by a bounded probability that the alignment is fortuitous.

Question of the proximity between protein sequences in the light of information theory
Since the optimized alignment score of two protein sequences allows an access to both the mutual information between proteins and an upper bound that the alignment is not fortuitous, one would expect that it is an accurate way to spatially organize proteins sets. A simple relation would be "the higher the mutual information, the nearest". There are three ways to assess the proximity between two objects a and b in a given space E [41]. The first is dissimilarity, a function f(a,b): a); the second is the distance per se, that is a dissimilarity such that the triangle inequality is respected: and the third is the similarity defined as a function f(a,b): a). Representing objects in a space is convenient using the notion of distance. When the optimal alignment is global, i.e. requiring that it extends from the beginning to the end of each sequence [32], it is theoretically possible to define a distance per se, that is to spatially organize the compared sequences [41]. However, from a biological point of view, global alignment algorithms are not reliable to assess homology of protein domains. Local align-ments are better suited, using scoring matrices to find the optimum local alignment and maximizing the sum of the scores of aligned residues [28,31]. In contrast with global alignments, local alignments do not allow any trivial definition of distances [41].
Although amino acid similarity is a function f(i, j): E × E → , owing to the local alignment optimization algorithms, the computed score is a function f(a,b): E × E → + , requiring the existence of at least one positive score in the similarity matrices. Thus, when constructing an alignment with the Smith and Waterman [33] method, the constraint that s(a,b)>0 (i.e. I(a;b) > 0) is imposed. This condition is consistent with proposition 1: if two sequences are homologous, knowledge about the first has to bring information about the second, that is to say, the mutual information between the two sequences cannot decrease below zero: I(a;b) > 0 (i.e. s(a,b) > 0). As a consequence, in the following geometric construction we sought a refined expression for the proximity of proteins.

Geometric construction of a configuration space of homologous proteins (CSHP) conserving mutual information
In a set of homologous proteins, any sequence a can be selected as a reference, noted a ref , in respect to which the others are compared. A geometric representation of objects relatively to a fixed frame is known as a configuration space (CS). In physics, a CS is a convenient way to represent systems of particles, defined by their positional vectors in some reference frame. Here, given n similar sequences, it is therefore possible to consider n references of the CSHP. In a given (CSHP, a ref ), each amino acid position aligned with a position in the a ref sequence, corresponds to a comparison dimension (CS dimension). Proteins are simply positioned by a vector, the coordinates of which are given by the scores of aligned amino acids. Gaps are additional dimensions of the CS. When considering that local algorithms identify the space of biological interest, i.e. a CSHP, the gap penalty is a parameter that maximizes the shared informative dimensions. Thus, given the amino acids mutual information, alignment optimization methods define the relative positions of proteins.  . This geometric stability is not observed with multiple alignments, which can be deeply modified by addition or removal of sequences.
In the CSHP, protein position is unaltered by additions or removals of other proteins. In practice, the construction of CSHP is therefore completely deduced from any all-by-all protein sequence comparison [45,46] and can be easily updated.

The q-dissimilarity, a proximity notion for a geometric representation of the CSHP
In the CSHP, the definition of a distance per se based on mutual information is reduced ad absurdum (For demonstration, see methods). To define a proximity function i) sharing properties of distance, i.e. increasing when objects are further apart, ii) deriving from similarity and iii) relying on mutual information, particularly the property Let s be a similarity, then q = e -s is a q-dissimilarity, named the 'canonical q-dissimilarity' associated to s. Accordingly, the TULIP theorem allows a statistical characterization of q(a,b) the canonical q-dissimilarity between two sequences a and b.

TULIP corollary 3
From the TULIP corollary 2, relation (14) is simply deduced: with Q(a,b*) being the random q-dissimilarity variable associated with S(a,b*). Given a (CSHP, a ref ), each sequence b aligned with a is characterized by a q-dissimilarity q(a,b). In geometric terms, b can be represented as a point contained in a hyper-sphere B of radius q(a,b).
The representation of a (CSHP, a ref ) shown in Figure 1 is therefore in conformity with all constraints listed earlier and can also serve as a Venn diagram for the setting of events realized following a continuous random variable Q(a,b*). When a is compared to itself, it is set on a hypersphere A of radius q(a,a), which is not reduced to one point. In the context of information theory, it is therefore possible to express that the proximity respects the property "q(a,a) ≠ q(b,b) is possible". Considering Figure 1, is the probability for a random sequence b* to be in the hyper-sphere B. In conclusion, the q-dissimilarity is therefore a proximity notion that allows a rigorous geometric description of the configuration space of homologous proteins, real or simulated, (CSHP, a ref , q).

Unification of pair-wise alignments theory, information theory, p-distance and q-dissimilarity in the CSHP model
A geometric space is a topological space when endowed with characterized paths that link its elements. Here, paths can be defined as the underlying evolutionary history separating sequences [11]. Given u the common unknown ancestor, then the divergence time t(a,b) is theoretically the summed elapsed times separating u to a and to b. Without any empirical knowledge of u, the simplest approximation for t(a,b) was sought as a function of the fraction of identical residues f id , thus of the p-distance. With the hypothesis of the molecular clock, this function can be given as equation (2), where the transmutation of a and b is a consequence of a Poisson process. By using relation (9) on the equivalence between score similarity and mutual information, then the fundamental postulate "the closer in the evolution, the more alike and conversely, the more alike, probably the closer in the evolution" can be reformulated:

Fundamental postulate
Given two homologous proteins a and b, the closer in the evolution, the greater the mutual information between a and b (i.e. the optimal computed score s(a,b)) and conversely, the greater the mutual information between a and b, probably the closer in the evolution.
Whereas the first part of the postulate is a consequence of the conservational pressure on mutual information, the second assertion founds the historical reconstruction underlying a set of biological sequences on statistical concepts. A corollary is that evolution of two homologous proteins is characterized by a loss of mutual information.
In the CSHP, this formulation of the fundamental postulate allows a novel mathematical formalization of the pdistance in probabilistic terms. Basically, the p-distance is the divergence observed between two sequences knowing that they share some features (the observed sequences a and b) and that they were identical before the speciation event (sequence u).
Looking back to equation (1), we can re-formulate f id in probabilistic terms, considering the fraction of shared features (identical sites) knowing the observed data and the existence of a common ancestor. Given two proteins a and b, let us consider the random variable Q(a,b*), defined in TULIP corollary 3. In (CSHP, a ref ), shown in Figure 1, one can define the probability law P{Q(a,b*)≤ρ} as the probability that the q-dissimilarity between b* and a ref is lower than ρ. The hyper-sphere of radius ρ contains therefore the b* random sequences sharing informative features with a accordingly. The probability p id/a that b* shares identity with a, knowing that the q-dissimilarity between b* and a is lower than that between the real sequences b and a, is: which is a probabilistic expression of f id in respect to the reference a ref . According to the Venn diagram in Figure 1: Using the Bayes theorem, equation (15) can be expressed as: In consequence: which can be expressed as

17
Assuming that substitution rates are independent of lineages [35], then random sequence models a* and b* are equivalent, that is to say Q(a,a*) ≈ Q(a,b*) and Thus p id/a , and symmetrically p id/b , provide a probabilistic expression of f id knowing the data, i.e. the observed mutual information between a and b expressed as Q(a,b).
Given two homologous sequences a and b, when their optimal score is s(a,b) ≥ µ + ψ with ψ being a critical threshold value depending on the score distribution law (See Methods for the demonstration for the critical threshold), owing to the TULIP corollary 2, we can state that p id/ a is bounded above: This expression can also be developed as: where µ 1 , σ 1 , µ 2 and σ 2 are the mean and the standard deviation of S(a,b*) and S(a,a*) respectively. The right term in relation (21) exhibits analogies with f id given by equation (1), showing that the pragmatic approach by Feng and Doolittle [19] could be supported and generalized in a theoretical elaboration.
Using the Poisson correction, an expression of t(a,b) is given as the linear combination of the two corrections of the p-distance deduced from p id/a and p id/b : with a* and b* the random variables corresponding to the shuffled sequences of a and b respectively. The sum of the logarithms corresponds to the product of the two probabilities, an expression of the hypothesis of independence of lineage. Interestingly, equation (22) provides an expression of the symmetric effect of time on the variations that independently affected a and b.
From relation (20), t(a,b) appears as a function of Z-score ratios. For any set of homologous proteins, it is therefore possible to measure a table of pair-wise divergence times and build phylogenetic trees using distance methods.

Reconstruction of protein phylogeny: first example, case study of the glucose-6-phosphate isomerase phylogeny
We compared the trees we obtained, called TULIP trees, to phylogenetic trees built using classical methods, for instance the popular PHYLIP [47] or PUZZLE-based [48] methods, termed here MAB trees (for multiple alignmentbased trees). Firstly, because MAB trees are constructed from multiple alignments, removals or additions of proteins modify the multiple alignments. Inclusion of sequences is considered as a way to improve the quality of multiple alignments and to increase the sensitivity of the comparison of distant sequences [49,50]. By contrast, the protein space used to build TULIP trees is not reordered when data sets are incremented or decremented (drawing of the TULIP tree may apparently change due to the tree graphic representation methods; nevertheless the absolute tree topology is not reordered). This remarkable property is due to both the geometrical construction by pairwise comparison and the convergence of the distance matrix elements estimated by equation (21). Indeed, the estimate of the right-hand term of equation (21) relies on a Monte Carlo method, after randomization of the biological sequences [39,44,51] and is therefore dependent on the sequence randomization model [52] and convergent in respect to the weak law of large numbers [53]. Convergence is proportional to , where numb rand is the number of randomizations. In the case studies presented here, we set num rand = 2000 (see Methods). By contrast, stability of MAB trees is sought by bootstrapping approaches and consensus tree reconstruction. MAB trees appear as the result of a complex learning process including possible re-adjustment of the multiple alignments after eye inspection pragmatically applied to assist the reconstruction. Alternatively, Bayesian analyses have been recently proposed for phylogenetic inference [54], estimating posterior probability of each clades to assess most likely trees. Still, in a recent comparative study, Suzuki et al. [55] and Simmons et al. [56] provided evidence supporting the use of relatively conservative bootstrap and jacknife approaches rather than the more extreme overestimates provided by the Markov Chain Monte Carlobased Bayesian methods. In the absence of any decisive methods to assess the validity of the trees obtained after such different approaches, no absolute comparison with the TULIP classification trees can be rigorously provided.
Whenever a TULIP classification was achieved on a dataset that led to a consensual MAB tree, both were always consistent. For example, Figure 2 shows the phylogenetic PHYLIP [47] and TULIP trees obtained for glucose-6- 1/ numb rand phosphate isomerases (G6PI). Phylogeny of the G6PI enzyme has been studied by Huang et al. [57] in order to demonstrate the horizontal transfer of this enzyme in the apicomplexan phylum due to a past endosymbiosis [57].
Owing to the neighbor-joining analysis used by Huang et al. [57] (see methods) Figure 2A shows that apicomplexan G6PI is "plant-like". The TULIP tree shown in Figure 2B is consistent with this conclusion. Interestingly, differences between the two trees are found only when the bootstrap values on the MAB tree are not strong enough to unambiguously assess branching topology.

Reconstruction of protein phylogeny: second example, case study of the enolase phylogenic incongruence
TULIP classification tree further helps in solving apparent conflicting results obtained with MAB methods. In a comprehensive study from Keeling and Palmer [36] the PUZ-ZLE-based reconstruction of the enolase phylogeny led to incongruent conclusions. Enolase proteins from a wide spectrum of organisms were examined to understand the evolutionary scenario that might explain that enolases from land plants and alveolates shared two short insertions. Alveolates comprise apicomplexan parasites, known to contain typical plant features as mentioned above, particularly a plastid relic. In this context, the shared insertion in apicomplexan and plant enolases (Figure 3) has been interpreted as a possible signature for Glucose-6-phosphate isomerase phylogeny some evolutionary relationship between apicomplexans and plants [58,59] and a likely sign of a lateral transfer. From the distribution of this insertion in enolases from several key eukaryotic groups, Keeling and Palmer [36] postulated that lateral transfer had been an important force in the evolution of eukaryotic enolases, being responsible for their origin in cryptomonads, Chlorarachnion and Arabidopsis. However, they could not conclude about alveolates, finding a conflict between the distribution of the insertion and the MAB phylogenetic position ( Figure 4A). The authors had to admit that lateral gene transfers failed to explain apicomplexa enolases, and were compelled to suppose that the lack of congruence between insertion and phylogeny could be because of a parallel loss of insertions in lineages, or to more complex transfers of gene portions.
Based on our theoretical model, we constructed the corresponding TULIP tree. TULIP trees given with BLOSUM 62 or PAM 250 matrices, Fitch-Margoliash or neighbor-joining methods led indistinctly to a unique tree topology ( Figure 4B). Separation of great phyla (Archaebacteria, Eubacteria, Diplomonads, Trypanasomes, Animals, Fungi and Amoeba) is recovered. A plant-like cluster is additionally reconstructed, in which a distinct separation occurred between {Rhodophytes ; Cryptomonads} and {Land Enolase phylogenic incongruence Figure 3 Enolase phylogenic incongruence. When aligned, the enolase region corresponding to amino acids 73-118 of the Oryza sativa gene, exhibit two insertions (red boxes) that are only present in land plants, charophytes and alveolates. In alveolates, these insertions are consistent with a horizontal gene transfer. However, to date, evolutionary reconstructions based on enolase sequences did not allow any phylogenetic branch gathering for these clades [36].
This topology corresponds to the observed distribution of the enolase short insertions and provides therefore a solution to the apparent enolase phylogeny incongruence: the phylogenetic position of alveolates is not in conflict with the distribution of enolase insertion and the apicomplexa enolase is possibly a consequence of a lateral transfer, like in cryptomonads.

Large scale phylogeny based on a CSHP built from massive genomic pair-wise comparisons
A CSHP containing large sets of protein sequences can be built after any all-by-all massive comparison providing Zscore statistics. Because the space elaboration is explicit, then quality of the mutual information conservation depends on the choice of the scoring matrix, the geometric positioning depends on the local alignment method, the homology assessment depends on the alignment score and probabilistic cutoffs and the phylogenetic topology on the choice of the stochastic law correction. Eventually, replicates. To build TULIP trees, Z-scores were estimated with 2000 sequence shuffling. Clades that contain a unique insertional signature (Land plants, green box; Charophytes and Chlorarachnion, blue box; Alveolates; yellow box) are not gathered in the MAB tree, as previously reported [36]. By contrast, in the TULIP tree, the phylogeny of enolase proteins is reconciled with the insertional signature detection in Land plants, Charophytes, Chlorarachnion and Alveolates.  genome-scale pair-wise comparisons [39,36]

Demonstration that distance of a protein to itself cannot be defined in the CSHP
In the simplest case, building a distance between amino acids (that would lead to distance between sequences) on the basis of computed similarity values would have to respect the condition: for i and j, two given amino acids and d the distance function. Using this condition in the proposition, any organization of the CSHP with a geometric distance is reduced ad absurdum.

Proposition
Building a distance between amino acids derived from the composed function d(i,j) = (φ ❍ s)(i,j), where s is a similarity function and φ a bijection, is impossible without a loss of mutual information. Moreover, two proteins from distinct organisms can have the same configuration, being like "twins", and d(i,j) = 0 does not imply i = j.

Proof
Condition (a) implies that φ (s(i,i)) = φ (s(j,j)) = 0. This equality imposes that s(i,i) = s (j,j) and, following equation (7) of main text, that I(i;i) = I(j;j). Considering for example tryptophan (W) and glutamic acid (E), if W occurs in a sequence, the mutual information gained about the occurrence of W at the aligned position would be the same as that gained in the case of E about the occurrence of E at the aligned position in the homologous protein. This statement is easily rejected on the basis of biochemical concerns. On one hand, aspartic acid (D) shares common biochemical properties with E, particularly a carboxylic acid, and easily substitutes in homologous sequences. By contrast W, exhibiting a unique biochemical feature, is less substitutable without altering the function. Thus the mutual information I(E;E) is necessarily lower than I(W;W). This that can be checked in scoring matrices such as BLOSUM 62 [30] where I(E;E) = 5, I(D;D) = 6 and I(W;W) = 11. Condition d(i,i) = 0 leads to an obvious loss of information. The second assertion of the proposition is obvious.

Determination of the threshold value ψ, for topological reconstructions in the CHSP based on pair-wise alignment score probabilities
An important basis of the reconstruction of a probabilistic evolutionary topology in the CSHP is based on the dem-onstration that, given S the random variable corresponding to the alignment scores of pairs of shuffled sequences and µ the mean of S, given two homologous sequences a and b, when their optimal score is s(a,b) ≥ µ + ψ (with ψ a critical threshold value depending on the score distribution law), owing to the TULIP corollary 2, we can state that p id/a is bounded above To the purpose of this demonstration, we considered the cumulative distribution function F(s) = P(S ≤ s), its derivative f(s) known as the probability density function defined as dF(s) = f(s)ds, and the positive delta function δ (s) = (s -µ) 2 (1 -F(s)). Since δ (s) = (s -µ) 2 (1 -F(s)) is null for s = µ and , the Rolle's theorem implies that ∃s 0 ∈]µ,+∞[ such as [60]; s 0 corresponds to a maximum of δ (s) and is therefore the solution of the one can express as The term corresponds to a continuous function. Interestingly, is known as the hazard function [61], that is the probability of s, per score unit (i.e. mutual information), conditional to the fact that the pair-wise alignment score is at least equal to s. The hazard function is also defined by . A critical hypothesis is that φ (x) function is strictly increasing and conversely that is strictly decreasing. Considering , equation 2 has only one solution s 0 and this solution is bounded above: In consequence, δ (s) reaches its maximum for a s 0 (s 0 ≤ µ + ψ) and it is strictly decreasing on ]µ + ψ ,+∞[. The estimation of s 0 is not trivial because it depends on the knowledge of the cumulative distribution function. Extensive studies provided experimental and theoretical supports for an extreme value distribution of alignment scores [31,43,44]. Using the extreme value distribution of type I, i.e. the Gumbel distribution [62], the cumulative distribution is given by with θ and β (β > 0) the location and scale parameters.
The probability density function g(s) is defined by dG(s) = g(s)ds. We observe with that . Using the Taylor's polynomial formula, i.e.

:
In consequence, for a Gumbel score probability distribution: A graphical determination of ψ from a Gumbel distribution is illustrated in Figure 5.
If a pair-wise alignment score of two sequences a and b is relatively high, that is s(a,b) ≥ µ + ψ, then the trivial inequality s(a,a) ≥ s(a,b) implies (s(a,b) -µ) 2 (1 - F(s(a,b))) ≥ (s(a,a) -µ) 2 (1 - F(s(a,a))) (h) that is to say From inequality (i), we deduce that p id/a is bounded above.

Construction of PHYLIP multiple alignment based trees and pair-wise alignment based TULIP trees
To build PHYLIP trees, multiple sequence alignments were created with ClustalW [63]. PHYLIP trees where constructed using the protpars and neighbor modules from the PHYLIP package [47] and the BLOSUM 62 substitution matrix. Bootstrap support was estimated using 1000 replicates. To build TULIP trees, for each couple of sequences a and b, alignment was achieved with the Smith-Waterman method and the BLOSUM 62 scoring matrices, using the BIOFACET package from Gene-IT, France [64]. We computed estimated z-scores z(a,b*), z(a,a*), z(a*,b), z(b*,b*), with 2000 sequence shuffling. For all computations, an estimation of the Gumbel parameters θ and β was made using the computed µ and σ of any S(a,b*) and the formula and θ = µ -βΓ'(1), where Γ'(1) ≈ 0.577216 is the Euler constant. In all computations, both Gumbel parameters were very close (in the case of enolases, mean(θ) = 35.04, SD(θ) = 0.12, mean(β) = 3.92, SD(β) = 0.08). As a consequence, the assumption Q(a,a*) ≈ Q(a,b*) was verified for any pair of sequences. We used the parameters to estimate µ = θ + βΓ '(1) (in the case of enolases, µ = 37.33), and µ + ψ ≈ µ + 10.5178 ≈ 47.85. As any pairs of computed scores are higher than this critical threshold, we used relation [20]. Estimation of evolutionary time was achieved according to equations [20] and [22]. Trees were constructed using Fitch-Margoliash and Neighbor-Joining methods [47].

List of abbreviations
CSHP, configuration space of homologous proteins, TULIP, theorem of the upper limit of a score probability


extfieldsofstudy: Biology Medicine

text:
Anaplasma phagocytophilum, Sardinia, Italy

To the Editor: Anaplasma phagocytophilum (formerly Ehrlichia phagocytophila), a tick-transmitted pathogen that infects several animal species, including humans (involved as accidental "dead-end" hosts), is the causative agent of human granulocytic anaplasmosis (HGA). It is a pathogen of veterinary importance responsible for tickborne fever of ruminants and for granulocytic anaplasmosis of horses and dogs (1,2). HGA was first described in the United States in 1994 (2) and is emerging in Europe (3). Although only 2 human cases have been reported in Italy (4), serologic and molecular findings have shown A. phagocytophilum infections in dogs and Ixodes ricinus ticks (5). Incidence, prevalence, and public impact of HGA and horse granulocytic anaplasmosis are, therefore, unknown for this geographic area. From 1992 to 1996, an average rate of 13.4 cases/year/100,000 inhabitants of tick bite–related fever of unknown etiology has been reported on the island of Sardinia, Italy, which is considerably higher than the corresponding national average value of 2.1 cases/year/100,000 inhabitants. Moreover, 117 cases of tick bite–related fever, whose etiology remains obscure, have been reported from 1995 to 2002 in the central west coast area of the island. Local newspapers occasionally report deaths as a result of tick bites, although no HGA-associated deaths have been documented in Europe.

This study investigated A. phagocytophilum in Sardinia. From 2002 to 2004, veterinarians based on the central west coast of the island were instructed to collect EDTA blood samples when a suspected case of tick bite–related fever was found at their clinics. A total of 70 blood samples were collected from 50 dogs and 20 horses that showed tick infestation and symptoms consistent with tickborne disease, such as fever, anorexia, jaundice (only in horses), anemia, myalgia, and reluctance to move. Genomic DNA was extracted from the buffy coat obtained by centrifugation of 2 to 4 mL of blood, as previously described (6). Furthermore, DNA was extracted from 50 Rhipicephalus sanguineus ticks removed from 30 dogs. Primers EphplgroEL(569)F (ATGGTATGCAGTTTGATCGC), EphplgroEL(1193)R (TCTACTCTGTCTTTGCGTTC), and EphgroEL(1142)R (TTGAGTACAGCAACACCACCGGAA) were designed and used in combination to generate a heminested polymerase chain reaction (PCR) for the selective amplification of 573 bp of the groEL gene of A. phagocytophilum. The final 50 μL PCR volume of the first PCR round contained 5 μL of the DNA extraction, primers EphplgroEL(569)F and EphplgroEL(1193)R, and HotMaster Taq DNA polymerase (5u/μL, Eppendorf) according to the manufacturer's basic protocol (Eppendorf AG, Hamburg, Germany). Heminested PCR was performed by using 5 μL of each of the first PCR products and primer EphgroEL(1142)R. To confirm the PCR diagnosis, amplicons were digested with the HindIII restriction endonuclease (predicted digestion pattern: 3 fragments of 525 bp, 21 bp, and 27 bp). Anaplasma phagocytophilum DNA was obtained from strain NCH-1 and used as positive control in PCR reactions. Sequences were obtained by cloning the PCR products into the pCR2.1-TOPO vector (Invitrogen S.R.L., Milan, Italy) and using the ABI PRISM Big Dye Terminator Cycle Sequencing Ready Reaction Kit (Applied Biosystems, Foster City, CA, USA), according to the protocols supplied by the manufacturers. Sequences ({"type":"entrez-nucleotide","attrs":{"text":"AY848751","term_id":"56463268","term_text":"AY848751"}}AY848751, {"type":"entrez-nucleotide","attrs":{"text":"AY848747","term_id":"56463260","term_text":"AY848747"}}AY848747) were aligned to the corresponding region of other species belonging to Rickettsiales by using ClustalX (7). Genetic distances among species were computed by the Kimura 2-parameters method by using MEGA, and were used to construct bootstrapped neighbor-joining trees (8).

Of 120 DNA samples, 1 tick, 3 dog, and 3 horse samples generated the predicted band of 573 bp representative of the groEL gene of A. phagocytophilum. HindIII digestions confirmed PCR diagnosis (Figure A1). Two different groEL sequence types were obtained from 1 dog and 1 horse and confirmed by BLAST (http://www.ncbi.nlm.nih.gov/Education/BLASTinfo/information3.html) queries as A. phagocytophilum groEL sequences (average identity 99%; average E value = 0), indicating that sequences did not reflect contamination. Bootstrapped neighbor-joining trees confirmed the identity of the new sequences obtained, which are closely related to HGA strains isolated in Europe and the United States (Figure).

The molecular approach applied in this study established A. phagocytophilum in an area of Sardinia characterized by a high prevalence of tick bite–related fever in humans and animal species. To our knowledge, this is the first evidence of A. phagocytophilum in Sardinian dogs and horses and the first documentation of infection in Italian horses caused by pathogenic strains. Therefore, these findings suggest the emergence of Anaplasma phagocytophilum in Italy. Ixodes ricinus ticks are indicated as vectors transmitting A. phagocytophilum in Europe. Although only 0.3% of 4,086 ticks collected in 72 sites of Sardinia (9) have been identified as Ixodes, other tick species are better represented on the island (Rhipicephalus, 67.2%; Haemaphysalis, 24.1%; Dermacentor, 4.9%). A. phagocytophilum in 1 Rhipicephalus sanguineus could indicate a role of this tick in the epidemiology of HGA. Finally, these data indicate the presence of a potential threat to human and animal health and suggest activation of further epidemiologic surveillance and controls.

Figure

Bootstrapped neighbor-joining tree of several species belonging to Rickettsiales and identification of the strains isolated during the study as Anaplasma phagocytophilum. Strains associated to Sardinian groEL variants are closely related to European and ...

Anaplasma phagocytophilum, Sardinia, Italy
To the Editor: Anaplasma phagocytophilum (formerly Ehrlichia phagocytophila), a tick-transmitted pathogen that infects several animal species, including humans (involved as accidental "dead-end" hosts), is the causative agent of human granulocytic anaplasmosis (HGA). It is a pathogen of veterinary importance responsible for tickborne fever of ruminants and for granulocytic anaplasmosis of horses and dogs (1,2). HGA was first described in the United States in 1994 (2) and is emerging in Europe (3). Although only 2 human cases have been reported in Italy (4), serologic and molecular findings have shown A. phagocytophilum infections in dogs and Ixodes ricinus ticks (5). Incidence, prevalence, and public impact of HGA and horse granulocytic anaplasmosis are, therefore, unknown for this geographic area. From 1992 to 1996, an average rate of 13.4 cases/year/ 100,000 inhabitants of tick bite-related fever of unknown etiology has been reported on the island of Sardinia, Italy, which is considerably higher than the corresponding national average value of 2.1 cases/year/ 100,000 inhabitants. Moreover, 117 cases of tick bite-related fever, whose etiology remains obscure, have been reported from 1995 to 2002 in the central west coast area of the island. Local newspapers occasionally report deaths as a result of tick bites, although no HGA-associated deaths have been documented in Europe.
This study investigated A. phagocytophilum in Sardinia. From 2002 to 2004, veterinarians based on the central west coast of the island were instructed to collect EDTA blood samples when a suspected case of tick bite-related fever was found at their clinics. A total of 70 blood samples were collected from 50 dogs and 20 horses that showed tick infestation and symptoms consistent with tickborne disease, such as fever, anorexia, jaundice (only in horses), anemia, myalgia, and reluctance to move. Genomic DNA was extracted from the buffy coat obtained by centrifugation of 2 to 4 mL of blood, as previously described (6). Furthermore, DNA was extracted from 50 Rhipicephalus sanguineus ticks removed from 30 dogs. Primers EphplgroEL(569)F (ATGGTATGCA-GTTTGATCGC), EphplgroEL (1193) R (TCTACTCTGTCTTTGCGTTC), and EphgroEL(1142)R (TTGAGTA-CAGCAACACCACCGGAA) were designed and used in combination to generate a heminested polymerase chain reaction (PCR) for the selective amplification of 573 bp of the groEL gene of A. phagocytophilum. The final 50 µL PCR volume of the first PCR round contained 5 µL of the DNA extraction, primers EphplgroEL (569)F and EphplgroEL(1193)R, and HotMaster Taq DNA polymerase (5u/µL, Eppendorf) according to the manufacturer's basic protocol (Eppendorf AG, Hamburg, Germany). Heminested PCR was performed by using 5 µL of each of the first PCR products and primer EphgroEL (1142)R. To confirm the PCR diagnosis, amplicons were digested with the HindIII restriction endonuclease (predicted digestion pattern: 3 fragments of 525 bp, 21 bp, and 27 bp). Anaplasma phagocytophilum DNA was obtained from strain NCH-1 and used as positive control in PCR reactions. Sequences were obtained by cloning the PCR products into the pCR2.1-TOPO vector (Invitrogen S.R.L., Milan, Italy) and using the ABI PRISM Big Dye Terminator Cycle Sequencing Ready Reaction Kit (Applied Biosystems, Foster City, CA, USA), according to the protocols supplied by the manufacturers. Sequences (AY848751, AY848747) were aligned to the corresponding region of other species belonging to Rickettsiales by using ClustalX (7). Genetic distances among species were computed by the Kimura 2-parameters method by using MEGA, and were used to construct bootstrapped neighbor-joining trees (8).
Of 120 DNA samples, 1 tick, 3 dog, and 3 horse samples generated the predicted band of 573 bp representative of the groEL gene of A. phagocytophilum. HindIII digestions confirmed PCR diagnosis (see Appendix Figure, available at http://www.cdc. gov/ncidod/eid/vol11no08/05-0085_ app.htm). Two different groEL sequence types were obtained from 1 dog and 1 horse and confirmed by BLAST (http://www. ncbi.nlm.nih. gov/Education/BLASTinfo/informa-tion3.html) queries as A. phagocytophilum groEL sequences (average identity 99%; average E value = 0), indicating that sequences did not reflect contamination. Bootstrapped neighbor-joining trees confirmed the identity of the new sequences obtained, which are closely related to HGA strains isolated in Europe and the United States (Figure).
The molecular approach applied in this study established A. phagocytophilum in an area of Sardinia characterized by a high prevalence of tick bite-related fever in humans and animal species. To our knowledge, this is the first evidence of A. phagocytophilum in Sardinian dogs and horses and the first documentation of infection in Italian horses caused by pathogenic strains. Therefore, these findings suggest the emergence of Anaplasma phagocytophilum in Italy.

Williamsia muralis Pulmonary Infection
To the Editor: Bacteria of the genus Williamsia are mycolic acid-containing actinomycetes of the suborder Corynebacterineae (1). This suborder also includes the genera Gordonia, Mycobacterium, Nocardia, Corynebacterium, Rhodococcus, Dietzia, Skermania, Tsukamurella, and Turicella (2,3). Within the genus Williamsia, only 2 species have been reported: Williamsia muralis, isolated from a daycare center (4), and W. maris, isolated from the Sea of Japan (5). One important aspect shared by both species is their apparent lack of pathogenicity, since they have been isolated only from environmental samples.
An 80-year-old woman, whose medical history included allergy to penicillin and high blood pressure, was admitted to the cardiothoracic intensive care unit at Juan Canalejo Hospital Complex in La Coruña, Spain, because of a loss of consciousness following an aortic valve replacement. Physical examination showed a systolic murmur and an echocardiogram showed aortic stenosis. Transaortic peak pressure was 100 mm Hg, and the aortic valvular area was 0.3 cm 2 . A biologic valve prosthesis (Mitroflow 21, Sorin Group Canada, Ltd., Burnaby, British Columbia, Canada) was inserted under the cardiopulmonary bypass.
Forty-eight hours later, the patient had paroxysmal atrial fibrillation and a temperature of 39°C, with severe hemodynamic and respiratory impairment. She was intubated and intravenous drugs were administered. Blood and urine cultures were requested. Central venous pressure lines were changed, and cultures were obtained. Empiric treatment with levofloxacin, amikacin, and teicoplanin was started for the patient. One of 2 blood cultures was positive for Staphylococcus epidermidis, as were cultures from femoral and jugular venous lines. Although considered a contaminant, we observed that S. epidermidis was susceptible to empiric antimicrobial drugs.
One week later, a chest radiograph showed bilateral alveolar infiltrates suggestive of pulmonary edema (Figure). To rule out infection, bronchoscopy and protected specimen brush were conducted. An unidentified gram-positive bacillus was cultured from the brush sample. Urine cultures were positive for Candida kefyr, but the patient showed no evidence of candidemia. An echocardiogram showed no evidence of infective endocarditis. Since the patient's condition did not improve, levofloxacin was replaced with imipenem, and treatment with fluconazole was initiated. However, the patient developed septic shock, adult respiratory distress syndrome, and oliguric acute renal failure, and died of multiple organ failure.
On direct examination, a Gram stain of the protected specimen brush sample showed numerous gram-positive bacilli. After incubation for 48 h in either an aerobic or capnophilic atmosphere, >1,000 CFU/mL were observed on Columbia agar plates containing 5% sheep blood (BD Stacker Plates, BBL, Franklin Lakes,


Sticking to the nuclear envelope

![][1]

LULL1 (red) draws torsinA (green) away from the nuclear envelope.

A mutation in torsinA, an AAA ATPase protein, causes DYT1 dystonia in humans. On [page 855][2], Goodchild and Dauer identify two protein binding partners of torsinA, one that localizes to the ER membrane and one

Sticking to the nuclear envelope mutation in torsinA, an AAA ATPase protein, causes DYT1 dystonia in humans. On page 855, Goodchild and Dauer identify two protein binding partners of torsinA, one that localizes to the ER membrane and one that spans the inner membrane of the nuclear envelope. The team hypothesizes that DYT1 is a nuclear envelope disease, and that identification of such protein complexes will provide a mechanistic probe into a poorly characterized region of the cell.
In wild-type cells, the majority of torsinA protein localizes to the lumen of the ER, while a small proportion associates with the nuclear envelope. In cells carrying the disease-associated mutation, the proportion of torsinA at the envelope increases. Furthermore, a torsinA mutant that cannot hydrolyze ATP, and therefore becomes trapped on its target protein, selectively localizes to the nuclear periphery.
To find out what torsinA binds at the nuclear envelope, the team screened a series of nuclear envelope proteins in cells expressing GFP-labeled torsinA. Overexpression of LAP1, a laminassociated protein, increased the proportion of perinuclear torsinA, and further analysis demonstrates an association between torsinA and the lumenal domain of LAP1. Based on sequence comparisons, Goodchild and Dauer identified a second torsinA binding protein, LULL1, an ER protein that has a single transmembrane domain and a long lumenal domain similar to that of LAP1. The torsinA-LULL1 association is also dependent on the lumenal region of LULL1.
The association between lamins, LAP1, and torsinA suggests that torsinA dysfunction may cause problems typical of A Exocytosis in small steps euroendocrine cells have two modes of exocytosis. In one mode, the vesicle fuses to the target membrane, and in the other-referred to as kiss-and-run-the vesicle touches the membrane and creates a small temporary pore. On page 929, Richards et al. show that the small synaptic vesicles in hippocampal neurons also use both mechanisms. The two modes of exocytosis differ in their kinetics and amplitude of the release, and may generate different postsynaptic responses.
To determine whether neurons have two pathways for exocytosis, the team loaded only a fraction of the synaptic vesicles with a fluorescent dye, a trick that allows them to resolve the activity of individual vesicles. They saw two types of exocytic events occurring: fast large releases of dye and slow small ones. The N laminopathies, including improper chromatin organization, transcriptional problems, or overall changes in nuclear architecture. Because the function of AAA ATPase proteins is knownthey typically use the energy of ATP hydrolysis to disassemble multiprotein complexes or unfold individual proteins-the researchers think they can use the complex as a launching point to learn about the function of the nuclear envelope.

LULL1 (red) draws torsinA (green) away from the nuclear envelope.
researchers conclude that the fast release events occur when a vesicle fuses with the plasma membrane, dumping its total dye content at a rate limited only by diffusion. On the other hand, the small release events appear to occur when a vesicle forms only a small pore in the plasma membrane. The speed of dye release in these small events is the same as that which occurs if the researchers use a toxin to poke a 1-2-nm hole in the membrane of a synthetic vesicle.
Richards et al. hypothesize that these different types of synaptic release may have a functional role in neural signaling. For example, it is known that when the vesicles of these hippocampal neurons fuse with the membrane and dump their whole load of neurotransmitter into the synaptic cleft, they excite the postsynaptic neuron. It could be that the dribbling release of neurotransmitter that occurs in the slow small exocytic events instead desensitizes the postsynaptic neuron. The data to confirm this are lacking, but several groups are now trying to figure out how to combine the electrophysiology and fluorescence experiments to find out.


extfieldsofstudy: Chemistry Medicine

s2fieldsofstudy: Biology Medicine

text:
Breakdown of accommodation in nerve: a possible role for persistent sodium current

Background Accommodation and breakdown of accommodation are important elements of information processing in nerve fibers, as they determine how nerve fibers react to natural slowly changing stimuli or electrical stimulation. The aim of the present study was to elucidate the biophysical mechanism of breakdown of accommodation, which at present is unknown. Results A model of a space-clamped motor nerve fiber was developed. It was found that this new model could reproduce breakdown of accommodation when it included a low-threshold, rapidly activating, persistent sodium current. However, the phenomenon was not reproduced when the persistent sodium current did not have fast activation kinetics or a low activation threshold. Conclusion The present modeling study suggests that persistent, low-threshold, rapidly activating sodium currents have a key role in breakdown of accommodation, and that breakdown of accommodation can be used as a tool for studying persistent sodium current under normal and pathological conditions.

Background
Accommodation is important for information processing in nerve fibers, as it determines whether, and how frequently, slowly-changing natural and artificial stimuli are translated into action potentials. Hill's theory of accommodation in nerve has been one of the most influential theories in this area [1]. A prediction of this theory is that a linearly rising current requires a certain critical slope in order to excite nerve fibers. Although this critical slope has been demonstrated in experimental preparations [2,3], it has not been found under normal physiological conditions [4,5]. Instead, nerve fibers have been shown to exhibit breakdown of accommodation; that is, a longduration slowly rising current excites nerve fibers at a nearly constant intensity no matter how slowly this intensity is approached [4,5]. A critical slope has only been found for depolarized nerve fibers, and Hill's theory of accommodation has been shown only to be applicable to such fibers [6]. Accommodation and breakdown of accommodation were the foci of several studies before the invention of the voltage-clamp, since prior to this innovation it was one of the few methods by which membrane kinetics could be studied. Since the invention of the voltage-clamp and later the patch-clamp some fifty years ago, the concept of breakdown of accommodation has been virtually absent from the scientific literature [7]. However, the biophysical mechanism responsible for breakdown of accommodation is still unknown; and as will be shown in this paper, a model that only contains transient sodium channels (i.e. currents that activate and deactivate rapidly in response to membrane depolarization) is unable to reproduce the phenomenon. Persistent (no inactivation) and late (slow inactivation) sodium channels have been identified in large dorsal ganglion neurons [8] and it has been found that these channels are needed for modeling latent addition in motor and sensory nerve fibers (i.e. threshold changes to short sub-threshold stimuli [9]). This suggests that persistent or late sodium channels are present in both motor and sensory myelinated nerve fibers and have fast activation kinetics that can initiate action potentials. The present study was undertaken to study the hypothesis that persistent sodium channels create a "threshold region" of membrane depolarization that cannot be exceeded without the generation of an action potential. Thus, it is suggested that persistent sodium channels are the cause of breakdown of accommodation.
The results in the present paper were based on a model of a space-clamped nerve fiber. This model included a persistent sodium channel based on the work of Bostock and Rothwell (1997) [9]. This channel was defined from the transient sodium channel with the following modifications: a) inactivation was removed (a persistent channel); b) the time-constant was slowed by a factor of two (timeconstant); and c) the kinetics was displaced so that the channel was activated at a membrane potential 20 mV more negative than is required to activate the transient channel (voltage shift) [9].

Model validation
The structure of the model and the choice of parameters allowed it to reproduce four sets of independent experimental data: threshold electrotonus, recovery cycle, latent addition and breakdown of accommodation (see Figure  1). The model was found to have a strength-duration time constant of 133.2 µs, which is similar to the experimental recorded value for the median nerve (139 ± 59 µs [10]). Furthermore, the model simulated breakdown of accommodation (see Figure 1D). The initial critical slope of the model was found to be 17.3 rheobase/s, which is lower than the experimentally recorded value (21.2 ± 3.72 rheobase/s) for ulnar nerves. The breakdown of accommodation seen in the model was likewise greater than observed in the ulnar nerve, and was closer to that observed in sensory nerve fibers [5]. The accommodation curve flattened out and remained near 2 rheobases when the time-constant of the current rise was greater than ~150 ms.

Breakdown of accommodation
In order to study the relationship between the properties of persistent sodium channels and breakdown of accommodation, one parameter at a time (number of channels, voltage shift, time constant) was changed and its influence on breakdown of accommodation was assessed (see Figure 2). When the voltage shift was decreased from -20 mV to -10 mV or -0 mV, it was still possible to create breakdown of accommodation by increasing the number of persistent sodium channels to 3.75% (-10 mV) or 15% (-0 mV) of the total number of sodium channels (see Figure 3A). However, for a voltage shift of -0 mV, the membrane potential did not return to the resting potential after the generation of an action potential (see Figure 3A).

Threshold responses to linearly rising stimuli
The threshold responses to linearly rising stimuli were significantly modified by the presence of persistent sodium channels (see Figure 4). (A threshold response is a response to a stimulus with intensity equal to the excitation threshold of the nerve fiber). Without persistent sodium channels, the threshold response to a linearly rising current of 20 ms duration did not occur at the end of the stimulus but had a latency of 5.04 ms (see Figure 4A). With 2.5% persistent sodium channels, the threshold response occurred at the end of the stimulus (see Figure  4B). This difference was found for all linearly rising currents tested that had stimulus durations in the range 1 ms to 200 ms; when the model had 2.5 % persistent sodium channels the threshold responses were always observed at the end of the stimulus, whereas without persistent sodium channels the longest latency of the threshold response was 5.04 ms (see Figure 4C). For both models, with and without persistent sodium channels, a non-linear response always occurred when the membrane was depolarized to a certain threshold value. This non-linear response initially occurred at the end of the stimulus, and with persistent sodium channels it resulted in an action potential. Without persistent sodium channels, it only resulted in an action potential when the stimulus intensity was sufficient for the response to occur with a latency of 5.04 ms or less.

Threshold electrotonus
A proportional relationship between the threshold change of a test stimulus and the underlying electrotonic changes in the membrane potential is a fundamental requirement for threshold electrotonus. Such a relationship was found for the model with 2.5% persistent sodium channels (see Figure 5). However, when the persistent sodium channels were removed from the model, the relationship between threshold and membrane potential broke down. This was tested with the conditioning current at an intensity of 40% for the model with persistent sodium currents, and 30.5% for the model without persistent sodium current. These two conditioning current intensities produced similar membrane depolarizations in the two models, which enabled the effect of the persistent sodium channels to be assessed. The relationship between threshold and membrane potential has also been found to break down when nerve fibers are depolarized because of a long-duration conditioning current or ischaemia [6]. Consequently, had the membrane depolarizations not been matched in the Comparison of the new model with experimental data for: A) threshold electrotonus [40], B) recovery cycle [41], C) latent addition [10], and D) accommodation curve [5] Figure 1 Comparison of the new model with experimental data for: A) threshold electrotonus [40], B) recovery cycle [41], C) latent addition [10], and D) accommodation curve [5]. In threshold electrotonus, a sub-threshold conditioning pulse of 100 ms duration is used to alter the threshold of a test stimulus delayed with respect to the onset of the conditioning pulse. In the recovery cycle, the nerve fiber is excited by a supra-threshold stimulus and the threshold of a test stimulus is determined at inter-stimulus intervals (T ISI ) of 2 ms to 100 ms. In latent addition, a short duration sub-threshold conditioning stimulus is used to alter the threshold of a test stimulus; the onset of the test stimulus is delayed with regard to the onset of the conditioning stimulus. In the accommodation curve, the threshold of stimuli of the form I S (1-e-tτ ) was determined, where τ was the time-constant of the current rise. In A, the bold line is the initial critical slope, which was estimated from the first four points in the accommodation curve where it is approximately a straight line. Experimental range: a) minimum and maximum of the experimental range, b) and c) mean ± standard deviation. D two models, the loss of the relationship between threshold and membrane potential may have been attributable to the greater membrane depolarization in the model without persistent sodium channels.

Discussion
We have used a model of a space-clamped motor nerve fiber to provide evidence for a link between persistent sodium currents and breakdown of accommodation. The model demonstrated that these channels might be the cause of breakdown of accommodation, as their inclusion enabled the model to reproduce the phenomenon (see Figure 1D). It also demonstrated that such channels are likely to be low-threshold and rapidly activating (see Figure 2). The low-threshold property is further supported by the fact that although breakdown of accommodation can be reproduced by high-threshold persistent sodium channels, in this case it results in an action potential that does not return to the resting potential (see Figure 3).   were therefore expected to modulate membrane excitability by amplifying and prolonging depolarization from a generator potential or an external electrode [8,11]. Indirect evidence has been obtained for the presence of such channels in both large diameter sensory nerve fibers and motor nerve fibers [9], and they can produce regenerative currents that facilitate action potential generation. Persistent sodium channels have been shown to amplify otherwise sub-threshold depolarization, thereby initiating action potentials [12]. Furthermore, acidification and alkalization within the physiological range have been found respectively to decrease and increase persistent and late sodium currents [13]. This pH-dependence of the late sodium current correlates well with experimental observations of breakdown of accommodation. Hence, breakdown of accommodation has been found to decrease during ischaemia and increase during hyperventilation [5]. Furthermore, when nerve fibers are depolarized with a polarizing current, there is a decrease in the threshold to triangular stimuli [14]. This suggests that it is not membrane depolarization per se that causes loss of breakdown of accommodation and the presence of a critical slope for slowly rising stimuli. In the present study, the loss of breakdown of accommodation is explained by loss of the persistent sodium current, such as would be caused by ischaemic depolarization due to acidification. Consequently, the present study predicts that the critical slope found by [2] was caused by ischaemic acidification and not membrane depolarization.

The effect of persistent sodium channels on threshold responses
In the present study, when the models with and without persistent sodium currents were stimulated by linearly rising currents, non-linear responses always resulted in action potentials when the model exhibited breakdown of accommodation. However, without breakdown of accommodation, non-linear responses only resulted in action potentials when they occurred within a "critical latency" from the onset of the stimulus. Without a "critical latency", which is the case with breakdown of accommodation, the threshold response occurred at the cessation of a long linearly rising stimulus. Consequently, the threshold for such stimuli is nearly constant regardless of their duration. However, when there is a "critical latency", the membrane potential needs to reach the voltage threshold within this "critical latency" for the nerve fiber to fire an action potential. The critical slope will then be proportional to the voltage threshold for which a non-linear response occurs divided by the "critical latency"; i.e. a decrease in "critical latency" results in an increased critical slope.

Model limitations
The present model simplifies existing knowledge of neuronal morphology and the distribution of ion channels. With regard to ion channels, only two potassium channels and two sodium channels were included in the model, but at least five distinct potassium channels [15] and three persistent and late sodium channels [11] have been identified, besides the classical transient sodium channel [16]. Unfortunately, current knowledge of the potassium and sodium channels in motor nerve fibers does not provide enough detail to allow modeling of them all. For example, the channel densities and kinetics are not known for all five potassium channels [15], and the kinetic data we have for the slow and fast potassium channels are likely to represent amalgamations of several channel species into single stereotypes [16]. Consequently, the present model is based on an amalgamation of distinct channels into stereotypes and the detailed geometrical structures of motor nerve fibers into a gross equivalent electrical circuit. The parameters of the model were based on experimental current-and voltage-clamp recordings whenever possible and the results obtained were found to be in line with experimental work. Consequently, these simplifications appear justified and therefore provide a basis for studying the biophysical properties of breakdown of accommoda-tion. This assumption is supported by previous work where models have provided insights into biophysical mechanisms [9,17,18]. However, the internodal leak resistance (R IL ) in particular was not based on experimental data, but was instead set by trial and error to a value that would enable the model to reproduce known experimental data. This approach was used since few experimental data on the internodal leak resistance are available. There are only modeling data on the periaxonal resistivity [19], but further modeling data suggest that the longitudinal conductance of the myelin sheaths has to be taken into account in determining the internodal leak resistance [20]. Consequently, an internodal leak resistance based solely on the width of the periaxonal space is likely to be an underestimate. The unknown resistivity of the periaxonal space presents further difficulties in obtaining a value for internodal leak resistance on the basis of experimental data alone. For these reasons we believe that the present approach was justified.

Alternative explanations of breakdown of accommodation
An alternative explanation for breakdown of accommodation could be the gating mode of the transient sodium channel [21]. The present paper follows the convention of assuming that activation and inactivation are two independent processes (i.e. the formalism of [22]). Today, it is known that activation and inactivation are inter-dependent, and that most transient sodium channels will go through an open state before entering an inactivated state [21]. This difference between the Hodgkin and Huxley formalism and recent knowledge of transient sodium channel function may have a synergistic role in breakdown of accommodation. Hence, a transient sodium channel with little inactivation before channel opening would not permit a critical slope and loss of breakdown of accommodation. However, this explanation remains unproven and would not change the conclusion of the present study, that persistent and late sodium channels can cause breakdown of accommodation. The interdependence of transient channel activation and inactivation may change the densities of persistent sodium channels needed for creating breakdown of accommodation, and thus there may be synergism between transient and persistent sodium channels.
A second explanation may be m-h overlap in the activation/inactivation kinetics of the transient sodium channel. For transient sodium channels, there is a region of membrane depolarization in which a persistent sodium current is generated l [23]. This is caused by channel activation while the membrane is still not sufficiently depolarized for all the channels to be inactivated, a phenomenon that has been termed m-h overlap. A theoretical study has demonstrated that the original squid axon model of Hodgkin and Huxley has breakdown of accommodation as a result of m-h overlap [23]. In this paper and other studies [9,24], persistent sodium channels are modeled as discrete channels. However, this does not imply that they are physically different from transient sodium channels. Three discrete persistent and late sodium currents have been identified on the basis of inactivation kinetics [11] in addition to the classical transient sodium current [16], but only one sodium channel Na v (1.6) has been found in the nodes of Ranvier in large peripheral nerve fibers [25]. This may suggest that persistent and late sodium currents are not generated specifically, but instead by transient sodium channels that operate in a gating mode with no or slowed channel inactivation. The modeling of persistent sodium current as created by persistent sodium channels does not provide evidence for the existence of such channels, only evidence that persistent sodium current can lead to breakdown of accommodation. Consequently, such persistent sodium current may be created by m-h overlap. However, in studies on persistent sodium currents, it has been argued that m-h overlap is not consistent with the observed kinetics [8,11]. Evidently, in mammalian nerve fibers, the persistent sodium current is most likely not generated by m-h overlap; but the study of [23] suggests that m-h overlap may be important for the persistent sodium current and breakdown of accommodation in squid axons.

Conclusion
The present modeling study has demonstrated that persistent sodium currents can create a "threshold region" for membrane depolarization that cannot be exceeded without the generation of an action potential. Thus, a persistent sodium current may be the underlying biophysical mechanism for the breakdown of accommodation to slowly rising currents, which are observed under normal physiological conditions in mammalian nerve fibers [4,5]. This suggests that accommodation curves can be used as a tool for studying persistent sodium currents under normal and pathological conditions.

Electrical model of a motor nerve fiber
The structure of the model of the space-clamped motor nerve fiber was based on previous models used for studying the accommodative properties of such fibers [9,26,27]. The present model represents a motor nerve fiber by the electrical equivalent circuit shown in Figure 6. The geometry of the node and internode was based on studies on the morphology of cat ventral spinal roots. The geometrical parameters were taken from cats of 1-11 years of age for a motor nerve fiber with a diameter of 14 µm (see Table 1). The nodal, internodal and myelin capacitances in the electrical equivalent circuit were calculated on the basis of these geometrical parameters and experimentally estimated capacitances per square micrometer (see Table 2 and Figure 6). The internodal leak resistance (R il ) and nodal resting potential were set by trial and error rather than calculated from geometrical and electrical parameters (see section entitled 'Validation', below).

Ionic currents
Five major ionic currents have been identified in myelinated nerve fibers as necessary for modeling a wide variety of experimental data: the transient sodium current (Na t ) for modeling the action potential [22], and the persistent sodium current (Na p ) for modeling latent addition [9] and the recovery cycle [24]. Fast (K f ) and slow (K s ) potassium currents have been shown to explain accommodation to depolarizing conditioning currents [28]. Accommodation to hyperpolarizing currents can be explained by a hyperpolarization-activated cation conductance (I H ), which is also thought to limit hyperpolarization in nerve fibers after they have conducted a train of impulses [28,29].
Transient and persistent sodium channels were included in the node, but following the work of [9] they were omitted from the internode for simplicity. The hyperpolarization-activated cation conductance was omitted from the model as it does not influence the response of nerve fibers to depolarizing stimuli [14]. Based on the work of [28,30,31]., the slow potassium current was included in the node as well as the internode. There is evidence for the localization of fast potassium channels in the paranode [32][33][34][35] As the paranode was not included in the present model, it was impossible to include fast potassium channels at this location. Instead, the approach used by [36] was applied and the fast potassium channels were included in the node. The ionic currents were described as being generated through membrane conductances (see Figure 6). The sodium con-   [50] ductances and slow potassium conductance in the node were based on single channel conductances and channel densities. Single channel conductances of 13pS and 8pS were used for the sodium channels and slow potassium channels, respectively [37]. The nodal densities for the sodium and slow potassium channels were set to 1000 channels/µm 2 [37] and 100 channels/µm 2 [38], respectively. The ion conductance of the fast potassium current was based on the work of [16], who found a fast potassium conductance of 15nS and a capacitive load of 1.4pF on the nodal membrane. The conductance of the fast potassium current was set from an estimate of the membrane area [16], which was based on the nodal capacitance in experimental data and the nodal capacitance per square micrometer [39].
The nodal resting potential was kept stable by a current leak to the internode, and the internodal resting potential was determined from this relationship. The internodal resting potential was kept stable by a small internodal sodium leak conductance. The nodal persistent sodium conductance was set by the fraction of nodal sodium channels that would be persistent. Therefore, the total number of nodal sodium channels was kept constant for all simulations.
Equivalent circuit for a space-clamped motor neuron Figure 6 Equivalent circuit for a space-clamped motor neuron. The model consisted of a node and an internode. Both the node and the internode contained non-linear current sources, which were calculated from equilibrium potentials and conductances. Channel types and maximum ionic conductances: node, transient sodium (Na t , 276nS), persistent sodium (Na p , 7.1nS), fast potassium (K f , 4.1nS), and slow potassium (K s , 17.4nS); internode, slow potassium (K s , 87.1nS) and leak conductance (L, 1.7nS). The linear parameters of the model were: C n , nodal capacity (0.22pF), C i , internodal capacity (379pF), C m , capacity of the myelin sheath (0.17pF), and R il , internodal leak resistance (41 MΩ).

Membrane kinetics
The non-linear membrane dynamics were based on human data [16]. The ionic current was given as: transient sodium current i Nat = G Nat m 3 h(E-E Na ), persistent sodium current i Nap = G Nap p 3 (E-E Na ), fast potassium current i Kf = G Kf n 4 (E-E K ), and slow potassium current i Ks = G Ks s(E-E K ).
The fractional activations (m, h, p, n and s) were given by the differential equation:  Table 3 for the constant: A, B, and C) and E is the membrane potential. The rate constants (α x and β x ) where scaled by appropriate Q 10 factors to a temperature of 37°C (see Table 3). All membrane kinetics were the same as those of [16] except the kinetics for the slow potassium current (see section entitled 'Validation', below). In order to allow the model to reproduce threshold electrotonus, it was necessary to modify the slow potassium channels. The kinetics were changed in order to slow the channel activation and to lower the fraction of open slow potassium channels at the resting potentials for the node and internode.

Validation
The model was validated with four sets of experimental data: threshold electrotonus [40], recovery cycle [41], latent addition [10] and accommodation curve [5] (slope and breakdown of accommodation). Threshold electrotonus is important as it provides insight into internodal conductances in human subjects in vivo, and it is promising for providing insight into disease mechanisms in neurological disorders [42]. In threshold electrotonus, subthreshold currents are used to alter the nodal and internodal membrane potentials. The change in threshold to a test stimulus is measured during the sub-threshold current, and this pattern of threshold alternations is termed threshold electrotonus [43]. The recovery cycle is a series of threshold fluctuations following an action potential. It is obtained by stimulating with a suprathreshold conditioning pulse and estimating the threshold with a subsequent test stimulus at various inter-stimulus intervals [43]. The threshold is usually tracked up to 200 ms after the conditioning pulse, during which time it goes through the absolute refractory period, relative refractory period, supernormal period and subnormal period. During the refractory and subnormal periods the threshold is increased, whereas it is decreased during the supernormal period [42,44]. Latent addition is obtained in the same manner as the recovery cycle [9,10]. The difference between the recovery cycle and latent addition is that the conditioning pulse is sub-threshold in latent addition but super-threshold in the recovery cycle. The strength-duration time constant τ was determined from the latent addition curve by fitting the function S2 = 100 -90e -s/τ to the simulated data, where S2 is the threshold of the test stimulus and s is the delay between the subthreshold conditioning stimulus and the test stimulus. Eleven delays, equally spaced between 0.0 ms and 1.0 ms, were used in this fit (see [10] for a more detailed description of the estimation of the strength-duration time constant using latent addition). The accommodation curve is a plot of the threshold current as a function of the timeconstant of current rise for exponentially rising stimuli [5]. Exponentially rising stimuli have the form I S (1-exp(t/τ)), where τ is the time-constant of current rise.
Five parameters were adjusted in order to fit the model to these experimental data: the nodal resting potential, the internodal leak resistance, the internodal slow potassium conductance, the nodal persistent sodium conductance and the kinetics of the slow potassium channel. Throughout the paper, modeling data are presented as superim- posed on the corresponding experimental ranges, a method taken from [24].

Implementation
The model was implemented in C and integrated by Euler's method with a time step of 2 µs. The model was interfaced with Matlab 6.0 as a mex function, and m-functions were written to estimate measurements of axonal excitability. The excitability measurements were based on a binary search algorithm, which determined the excitation threshold with an accuracy of 0.1pA. An action potential was identified if the nerve fiber was depolarized to -30 mV with a rate of rise of more than 60 mV/ms. Stimulation was achieved by an intracellularly-injected current in the node.