extfieldsofstudy:
Biology
Computer Science
Medicine

text:
Detection of prokaryotic promoters from the genomic distribution of hexanucleotide pairs

Background In bacteria, sigma factors and other transcriptional regulatory proteins recognize DNA patterns upstream of their target genes and interact with RNA polymerase to control transcription. As a consequence of evolution, DNA sequences recognized by transcription factors are thought to be enriched in intergenic regions (IRs) and depleted from coding regions of prokaryotic genomes. Results In this work, we report that genomic distribution of transcription factors binding sites is biased towards IRs, and that this bias is conserved amongst bacterial species. We further take advantage of this observation to develop an algorithm that can efficiently identify promoter boxes by a distribution-dependent approach rather than a direct sequence comparison approach. This strategy, which can easily be combined with other methodologies, allowed the identification of promoter sequences in ten species and can be used with any annotated bacterial genome, with results that rival with current methodologies. Experimental validations of predicted promoters also support our approach. Conclusion Considering that complete genomic sequences of over 1000 bacteria will soon be available and that little transcriptional information is available for most of them, our algorithm constitutes a promising tool for the prediction of promoter sequences. Importantly, our methodology could also be adapted to identify DNA sequences recognized by other regulatory proteins.

Background
Adaptation is essential to the survival of any biological organism and requires appropriate transcriptional regulation to modulate gene expression profiles. In prokaryotes, RNA polymerase (RNAP) is responsible for the transcrip-tion of all genes. However, promoter recognition is effected by an interchangeable sigma (σ) factor that associates to RNAP, directs the newly formed holoenzyme to a promoter and contributes to transcription initiation [1]. Transcription levels can be further modified by additional regulators (activators and repressors) that affect the recruitment or the activity of RNAP holoenzymes at various promoters [2]. A common feature of σ factors and transcriptional regulators is their ability to recognize specific DNA patterns in order to modulate gene expression. It is presumed that, as a result of evolutionary pressure, these regulatory sequences were selected upstream of some genes or operons and excluded from the rest of the genome.
Bacterial genomes usually encode many σ factors. Of these, the principal σ factor would be responsible for the expression of housekeeping function genes. The remaining σ factors are thought to direct the expression of genes required for specialized functions such as stress responses or sporulation [1]. σ factors can also be classified according to their structural homology to either σ 70 or σ 54 of Escherichia coli. σ 70 -related factors, which constitute the vast majority of known σ factors, are composed of two major DNA-binding domains capable of sensing a certain spacing range when associated to RNAP [3,4]. These σ factors usually recognize two DNA boxes (herein referred as the promoter) of approximately six base pairs (bp) located roughly at 10 and 35 bp upstream of the transcription start site (TSS). Spacing between these two boxes generally ranges from 16 to 20 bp [5]. σ factors similar to σ 54 also recognize two DNA boxes in the promoter region. However these elements are located approximately 12 and 24 bp upstream of the TSS. Other major differences between σ 54 and σ 70 family members are the ability of σ 54 to bind DNA in absence of RNAP and the requirement of an isomerization step by an activator to render σ 54 -containing holoenzymes processive.
σ factors can tolerate a variety of mismatches from their consensus sequence. For example, a typical E. coli σ 70 promoter sequence contains two mismatches within both the -35 and -10 hexanucleotide elements [6]. However, there is generally a direct relationship between promoter strength and the similarity to the corresponding consensus sequence [7]. Variations over three orders of magnitude have been reported in σ 70 -dependent promoter strength in E. coli [8]. In some cases, an extended -10 promoter box may be observed and may substitute for the absence of a clear -35 element. Extended -10 promoter boxes were reported to be present in 20% of promoter sequences in E. coli [6] and 45% in Bacillus subtilis [5].
A variety of techniques have been used to identify TSS and to characterize σ factor-DNA interactions. However, the formal identification of promoters by molecular methods can be tedious and is currently not amenable to genomewide applications. Consequently, it is important to develop algorithms that can rapidly and accurately evaluate the presence of promoters, without the need for exten-sive biochemical studies. Current algorithms for promoter detection, typically developed for a specific bacterium, exploit different characteristics of promoter sequences. Some approaches are based on sequence representation or statistical overrepresentation. Other methodologies have also been described for the detection of DNA motifs in sets of regulatory sequences [9][10][11][12] or by comparing the upstream regions of orthologous genes from different species [13][14][15][16][17]. A method based on the weaker stability of the DNA double-helix in promoter regions was also recently used to identify promoter regions [18]. However, most of these procedures are not suitable for the identification of precise prokaryotic promoters because of the inherent variability in promoter sequences and because they do not allow variable spacers between two DNA motifs.
Sequence representation strategies designed for promoter identification are usually based on a prior knowledge of some characterized sequences. These algorithms are thus trained to recognize sequences that are similar to a previously defined representation of a promoter. This approach was first used by Galas et al. [19]. As reported by Stormo, numerous false positives (FP) are however obtained with this strategy [20]. For example, allowing two mismatches in the σ 70 consensus -10 hexanucleotide produces roughly one hit per 30 nucleotides (nt) in the complete genome of E. coli. A more accurate representation of DNA-binding motifs consists of position-specific weight matrices (PSWM) [21], and online tools such as Virtual Footprint [22] are available to facilitate their analysis in the context of bacterial gene expression. Nonetheless, searching for full E. coli σ 70 consensus promoter sequences using more flexible mismatch restrictions offered by PSWM also generates a vast amount of hits [23]. More recently, Huerta and Collado-Vides used a PSWM-derived methodology and detected approximately 15 putative promoters/100 nt in IRs [24] By adding several constraints such as grouping sequences and filtering with the distance from the start codon, they achieved a sensitivity of 86% with an average of 1.88 putative promoters/100 nt. Several groups have also used general neural networks but no significant improvements have been achieved over the PSWM [25]. Hidden Markov Models (HMM) have also been trained to identify promoter sequences recognized by the principal σ factor in B. subtilis [26] and Campylobacter jejuni [27]. A learning approach based on a Support Vector Machine (SVM) employing a variant of the mismatch string kernel was also recently described [28]. Importantly, all above-mentioned approaches depend on a previously established or trained description of promoters, and were not designed to function with organisms for which promoter information is insufficient.
Statistical overrepresentation approaches can identify short DNA sequences that are present more frequently in a subset of sequences than what would be expected by chance according to the background distribution. Using such a procedure, Vanet et al. have proposed a description of the promoter sequences recognized by the principal σ factor of Helicobacter pylori from different sets of IRs [29]. More recently, the MITRA algorithm, which also evaluates the spacing between promoter boxes and the positional bias from the start codon, was applied to 20 bacterial genomes. Four of these genomes generated statistically strong signals possibly corresponding to principal σ factor-dependent promoter sequences, including the ones from H. pylori and B. subtilis [30]. Using a different approach, the principal σ factor consensus sequence was identified among over-represented motifs in B. subtilis, although the methodology was not designed especially for that purpose [31]. The latter study was based on the method of Li et al., designed to identify regulatory protein binding sites in E. coli [32]. It has been noticed that the E. coli σ 70 consensus sequence was not identified by this or other approaches, a failure that was attributed to the greater variability of promoter sequences within this organism [32]. A similar method was also unable to distinguish a motif related to the principal σ factor promoters in the complete genome of Streptomyces coelicolor [33]. In general, although some statistical approaches had limited success, these methods do not seem appropriate in their current form for the identification of promoter sequences in a variety of organisms.
In this paper, we describe a novel approach based on matrices representing the genomic distribution of hexanucleotide pairs, and designed to predict precise promoter sequences using any annotated prokaryotic genome. This approach can be applied to organisms for which almost no transcriptional data is available, without the need for extensive biochemical characterization. The strategy is based on the observation that, although promoter sequences can vary for every σ factor and according to the GC content of each genome [34], promoters are over-represented in IRs relative to the whole genome. Since this bias appears to be conserved throughout evolution, the characteristic distribution of promoter sequences is thus used to identify promoters in a variety of prokaryotic organisms. Briefly, a score is calculated based on the similarity between a matrix representing the genomic distribution of most promoter sequences reported in the literature and a matrix representing the genomic distribution of a putative promoter sequence. A Z-score is next calculated according to the background. To assess the validity of our method, over 680 characterized promoter sequences from ten genomes were gathered from databases and from the literature, and tested using various sta-tistical indicators. Experimental validations of promoter prediction also supported our approach.

Genomic distribution of regulatory sequences is biased and conserved
Transcription initiation is an important step in the regulation of most bacterial cell processes [35,36]. For this reason, it is thought that transcription factor binding sites have been evolutionary selected in some IRs, and generally excluded from the rest of prokaryotic genomes in order to avoid aberrant gene expression. Since regulatory sequences can often tolerate a variety of mismatches with respect to the corresponding DNA motif consensus sequence and still be efficiently bound by their cognate DNA binding protein, these degenerated sequences are also believed to be subjected to a similar selective pressure. Therefore, the genomic distribution of transcription factor consensus DNA binding sequences and close derivatives is expected to be biased towards IRs possibly involved in transcription initiation relative to the whole genome. To test this latter hypothesis, we have analyzed the genomic distribution of consensus sequences representing, respectively, the recognized sequences of a principal σ factor, an alternative σ factor and a transcriptional regulator of three bacterial species (Figure 1 and data not shown). The global distribution of these DNA sequences is reported in "distribution matrices", which consists of tables of dimension 4 × 4 where lines and columns correspond respectively to the number of mismatches in the proximal and distal boxes of a particular DNA sequence. For every matrix cells, the proportion of hits in putative promoter-containing IRs (P) relative to the whole genome (G) was calculated, with respect to a defined spacing range between the two hexanucleotide boxes forming the evaluated sequences. Light colored cells correspond to low P/G ratios while darker cells represent higher frequencies. For example, the genome-derived distribution matrix of the E. coli σ 70 consensus promoter sequence shows that 65% of genome sequences bearing one mismatch in the -35 element (distal box) but a perfect -10 (proximal) box are localized in P regions relative to the whole genome. At three mismatches in each box (50% degeneracy of the analyzed sequence), about 16% of the hits are found in P regions, which roughly corresponds to the proportion of IRs in the genome. Such sequences are thus considered randomly distributed in the genome. Results presented in Figure 1 thus demonstrate that transcription factor consensus DNA binding sites and related mismatch(es)-containing derivatives tend to be preferentially localized inside IRs. Moreover, the genomic distribution of these sequences shows related patterns in all tested bacterial species ( Figure 1).

Principal σ factor promoter sequences
Because principal σ factor promoter sequences are readily available and relatively abundant in many bacterial species, we decided to further explore their genomic distribution. Figure 2 (upper and middle rows) shows genomederived distribution matrices of experimentally identified promoter sequences from E. coli, B. subtilis and M. tuberculosis. As expected, these matrices contain patterns similar to those observed with the corresponding principal σ factor consensus sequences ( Figure 1). Importantly, the genomic distribution of sequences from the coding region of the rpoB genes from the same bacterial species did not show any enrichment inside of P regions ( Figure 2, last row and data not shown). These results suggest that transcription factor DNA binding sites from various bacterial species have a genomic distribution significantly different from that of non-regulatory sequences.

Synthetic matrices
From these observations, we hypothesized that a particular genome-derived distribution matrix may appropriately represent several promoters. Hence, it could be possible to calculate a score reflecting the similarity between the distribution matrix of a typical promoter and the genomederived distribution matrix of any hexanucleotide pair from the same organism. A high score would indicate a strong probability that the tested sequence is also a promoter. However, the best reference matrix should not necessarily be obtained from an existing promoter sequence. In fact, an interpolated matrix could indeed offer more flexibility and be much more effective. We therefore decided to synthetically generate distribution matrices according to the values observed in each cell of the genome-derived distribution matrices for all available experimentally identified principal σ factor-dependent Genome-derived distribution matrices generated for consensus sequences from transcription factors of three different organ-isms The name of the transcription factor is identified above each matrix. The mismatch number of each cell is indicated on both sides of matrices. The analyzed consensus sequence is shown under each matrix, along with the allowed spacing range. The first row corresponds to principal σ factors, the second row to alternative σ factors, and the third row to transcriptional regulators.
promoter sequences (see Additional file 1: A schematic description of the procedures used in this work). A range of ratios was next determined for each cell, resulting in over 248 million different synthetic matrices (see Additional file 2: Detailed information on the generation of synthetic matrices). Synthetic matrices are thus not affiliated to any hexanucleotide pairs, but are rather produced from the genome-derived distribution matrices of experimentally identified promoters. Moreover, it could be possible to identify synthetic matrices suitable to detect promoter sequences for a specific organism ("organismspecialized matrix"), and perhaps for all bacteria ("general matrix").

Using distribution matrices to detect promoters
To measure the ability of synthetic matrices to correctly identify promoter sequences, we tested them on 684 characterized principal σ factor dependent promoter sequences from ten bacterial species (Table 1). Complete IRs containing these promoters were extracted from their respective genomes along with 30 nt on each side. These "enlarged IRs" were next scanned using two hexanucleotide windows with respect to the allowed spacing ranges. For each hexanucleotide pair, a genome-derived distribution matrix was generated and used to calculate a score based on its similarity to the synthetic matrix under evaluation (see Additional file 3: Detailed example of score Specific examples of genome-derived distribution matrices generated for characterized principal σ factor-dependent promoter sequences and non-regulatory sequences from three different organisms Figure 2 Specific examples of genome-derived distribution matrices generated for characterized principal σ factordependent promoter sequences and non-regulatory sequences from three different organisms. (A) E. coli, (B) B. subtilis, and (C) M. tuberculosis. Gene names and experimentally identified promoter sequences, as reported in the literature, are indicated above each matrix. The first row corresponds to characterized promoter sequences closely resembling to the proposed consensus. The second row presents experimentally identified promoters containing more mismatches relative to the proposed consensus. The last row shows distribution matrices of hexanucleotide pairs with approximately 3 mismatches position per box, which were extracted from the middle of the rpoB coding sequence (bona fide non-promoter sequences). calculation). Different metrics have been tested and the selected one provided the best outcome (data not shown). Since the -10 box has a more important role in transcription [37], the mean of scores sharing the same -10 box was calculated. This latter averaged score had to be greater than the threshold to consider the corresponding sequence as a candidate promoter. Furthermore, since the distance between the TSS and the proximal hexanucleotide box can slightly vary (which may cause the experimental identification of -10 boxes to be inaccurate), signals located at ± 4 nt of an experimentally identified proximal hexanucleotide box were considered as true positives (TP, see Methods). Indicators such as sensitivity, specificity, precision, and performance were calculated for each synthetic matrix [38]. We chose to present the FP rate per 100 nt rather than the precision indicator to facilitate comparison since the precision indicator does not take into account the length of analyzed regions (Table 1). In contrast with other studies [24,39], the analyzed regions were not limited to a maximum length, longer regions producing more FPs.
"Organism-specialized" synthetic matrices (selected on the basis of the performance and sensitivity indicators) gave interesting results for each tested organism (Table 1). For instance, the synthetic matrix #113362653 identified almost 60% of promoters among the set of 148 characterized promoter sequences from B. subtilis with approximately one FP/100 nt. Overall, the sensitivity of the best matrix for each organism ranges from 29.4% to 90.9% with 0.53 to 1.42 FP/100 nt (Table 1). Amongst the FP, some could be uncharacterized real promoter sequences. Performance, precision and specificity indicators ranged respectively between 4.6-23.8%, 5.2%-24.4% and 98.6-99.5% (data not shown). Cross-validation tests were also conducted with E. coli and B. subtilis promoter datasets and matrices very similar to the organism-specialized syn-thetic matrices were identified. Moreover, the sensitivity and FP rate of these cross-validation matrices were comparable to the organism-specialized matrices obtained using complete datasets, demonstrating the robustness of the approach and suggesting that the various organism-specialized matrices are appropriate (see Additional file 4: Three fold cross-validation results).
Interestingly, several synthetic matrices identified a significant fraction of promoter sequences in all tested genomes. We thus undertook to find a matrix that would be the best compromise for promoter detection in all bacteria. This was achieved by summing up the relative performance indicator of each matrix on all genomes. The best "general" synthetic matrix (#45012859) ( Figure 3A and Additional file 3) usually showed a slight decrease in sensitivity relative to organism-specialized synthetic matrices (Table 1). Nonetheless, the sensitivity of this general synthetic matrix ranged from 14.7% to 90.9% with 0.70 to 1.36 FP/100 nt (Table 1), which is comparable to other approaches for an equivalent FP rate (see Additional file 5: Comparison with other bacterial promoter prediction approaches). Furthermore, we have replaced our scoring function by a function based on PSWM. The results obtained with E. coli and B. subtilis normalized PSWM promoter datasets, showed a significantly decreased sensitivity and a lower FP rate in our design, suggesting that the genomic distribution of promoter sequences is more flexible, and hence more portable between species than current promoter sequence models (see Additional file 6: Comparison with normalized PSWM scoring function).
The scores of all possible hexanucleotide pairs for a specific enlarged IR can be presented in a graph. As an example, Figure 3 shows graphs obtained using the enlarged IRs containing the characterized promoter sequences pre- For each organism, the GC content, the number of gathered characterized promoter sequences, the number of distinct IRs containing these promoters, the total number of nt in these IRs, and the allowed spacing between hexanucleotide pairs are indicated. The sensitivity and FP rate obtained with the best synthetic matrix (organism-specialized) and the general matrix are also shown.
sented in Figure 2 for E. coli ( Figure 3B), B. subtilis ( Figure  3C) and M. tuberculosis ( Figure 3D). The rpoB scanning is also presented for E. coli since identical results were obtained using the B. subtilis and M. tuberculosis sequences (Figure 3 and data not shown). Interestingly, Figure 3C (Bsu-lonA) shows one of many examples where a -10 promoter box (consensus TATAAT) is coupled to a ribosome binding site (RBS, consensus AGGAGG).
To assess the significance of our results, the analyses were repeated on two types of shuffled genome sequences. While retaining all the original genome information, the shuffling destroys its structure, which is used by our methodology to identify promoter sequences. The first shuffling procedure was accomplished by repositioning mononucleotides one region (gene or IR) at a time, thus keeping the AT bias intact in each IR [40]. The second shuffling was performed independently of gene annotations, thus dispersing the GC content uniformly through a genome. We surmised that most regulatory sequences, and by extension their genomic distribution, would be affected differently by these procedures. The overall sensitivity obtained with shuffled genomes should thus be decreased. Indeed, genome-derived distribution matrices were completely different if calculated from the real genome or from shuffled genomes (data not shown). As expected, the second shuffling was much more detrimental to the observed sensitivity. For instance, the calculated sensitivity with the general synthetic matrix dropped respectively from 31% (intact genome) to 14% and 0.6%  (shuffled genomes) for E. coli, and from 50% to, respectively 17% and 0.5% for B. subtilis. The same trend has been observed for other genomes (data not shown). Shuffling performed with longer nucleotides blocks than mononucleotide gave intermediate results (data not shown).

Coupling sequence information to genome distribution
For some bacteria, consensus promoter sequences representing the sequences recognized by principal σ factors are well characterized. For example, the consensus promoter sequence for B. subtilis and E. coli is TTGACA < >TATAAT [5]. We have thus examined the possibility of decreasing the occurrence of FPs by introducing a very simple sequence-dependent filter. Approximately half FPs were eliminated by allowing up to three mismatches in the -35 (distal) box and up to two mismatches in the -10 (proximal) box (the latter being more conserved in these organisms). For instance, a sensitivity of 49% with 0.50 FP/100 nt for the specialized matrix and 44% with 0.44 FP/100 nt for the general synthetic matrix was obtained for B. subtilis. Since fewer FPs were observed, the threshold could be lowered, which can improve the sensitivity of the algorithm. For example, the reduction of the IR specific threshold (see Methods) from three to two standard deviations upon sequence filtering gave a sensitivity of 58% with 0.80 FP/100 nt in B. subtilis using the general synthetic matrix. Figure 4 illustrates how the application of a filter may allow the identification of an otherwise missed promoter.

Prediction and biochemical validation of promoters
In addition to the IRs used to evaluate synthetic matrices, we have applied our promoter detection algorithm on the remaining IRs from the ten genomes mentioned in Table  1. Some predicted promoters from B. subtilis, E. coli and M. tuberculosis were chosen for further investigation using a biochemical approach. RNA was extracted from exponentially growing cultures and TSS were determined by primer extension analyses. The seven candidate promoters for which a transcription signal was observed matched the predicted sequence ( Table 2). The distance between the observed TSS and the potential -10 box was in all cases either 5 or 6 nt. These results demonstrate that our algorithm represent a valuable tool for the identification of prokaryotic promoter sequences. Complete prediction lists for these three organisms, and others, are available upon request.

Discussion
In this work, we present evidences suggesting that regulatory sequences and their close derivatives have a biased distribution pattern for IRs that may support transcription initiation. Furthermore, our data support the idea that the preferential location of regulatory sequences is shared between bacterial species. In order to clearly demonstrate the potential of genomic distribution as an indicator of DNA motif function, we have developed an algorithm that can identify a significant fraction of principal σ factor-dependent promoters in any prokaryotic organism, using only a genome annotation and a synthetic matrix (i.e. the general matrix, which was obtained from a training set composed of experimentally identified promoters from ten bacterial species). Promoter predictions were also made and experimentally verified, thus highlighting the potential of our approach for promoter identification in various prokaryotic organisms. Overall, our strategy yielded results similar to those from other studies considering an equivalent amount of FP/100 nt (see Additional file 5: Comparison with other bacterial promoter prediction approaches). However, our algorithm took advantage of a yet unexploited concept, can be used in a wide variety of organisms, required almost no previous knowledge of promoter sequences to be effective, and can be combined with other methodologies. The fact that our general matrix allowed the detection of more promoters in all tested genomes relative to the E. coli and B. subtilis principal σ factor normalized PSWM (see Additional file 6: Comparison with normalized PSWM scoring function) also supports the idea that genomic distribution of pro- The application of a simple sequence-dependent filter and the subsequent diminution of the threshold allows the detection of an otherwise missed promoter Figure 4 The application of a simple sequence-dependent filter and the subsequent diminution of the threshold allows the detection of an otherwise missed promoter. Graphs reporting scores calculated for all hexanucleotide pairs located in the enlarged IR upstream of the B. subtilis ftsA gene. (A) No filter applied. (B) Sequence-dependent filtering (3 mismatches and 2 mismatches allowed respectively in the -35 and -10 boxes relative to the B. subtilis principal σ factor consensus sequence) followed by the reduction of the region-specific threshold from 3 to 2 standard deviations above of the mean.
moter sequences is easier to transfer across organisms with regards to current promoter sequence models.
Although our approach is based on the genomic distribution of hexanucleotide pairs rather than a direct sequence evaluation, it is still important to know the approximate spacing range that is tolerated by a σ factor to efficiently detect the corresponding promoter boxes (data not shown). However, this distance appears to be very similar in most bacteria. The spacing range limitation restrains putative promoter signal contamination by irrelevant hexanucleotide pairs.
An important assumption in our method is that all promoter sequences share a related genomic distribution pattern. However, it is possible that some promoters fall in distinct biological categories or slightly differ between bacterial species. As a consequence, specific matrices could be more adapted to different promoter types. For example, a matrix could be more suitable for promoters containing an extended -10 promoter box. Similarly, very weak promoter sequences could bear an altered distribution pattern when compared to strong promoters.
Another important consideration in our study is the relatively small size of prokaryotic genomes. Since many of these contain only a few million bp, some hexanucleotide pairs are particularly absent from IRs and/or from the whole genome. Therefore, blanks (or 0) in mismatch containing cells can be found in some genome-derived distribution matrices, thus strongly altering the resemblance with the synthetic matrix. Nonetheless, a few promoter sequences containing a blank cell are identified by our approach (see Mtb-rrs in Figures 2C and 3C), although most of them are not (data not shown). Since the score is calculated from the mean of hexanucleotide pairs sharing the same proximal box, a blank occurring only at a particular spacing may not be too detrimental to the overall score of a sequence.
An additional possibility to explain our inability to identify some promoters is that, although present in the genome, some mismatch combinations may be relatively rare, and a small variation in absolute numbers may have a significant impact on P/G ratios. In addition, some promoter sequences contain nt triplets corresponding to codons frequently used in translation, which may flatten their distribution bias for IRs. For instance, 19% of the TP and 44% of the FN hexanucleotide pairs of E. coli evaluated with the organism-specialized matrix (22% and 41% respectively with the general synthetic matrix) include a triplet which is used as a codon more frequently than average. Thus, almost half of the FN in E. coli seem to have an impaired distribution profile because of the inclusion of at least one frequent codon.
In spite of the fact that our algorithm was designed for fully sequenced and annotated genomes, preliminary tests suggest that a genomic distribution calculated from a closely related organism can be used as a reference with similar results (data not shown). Similarly, errors in genomes annotation could theoretically have an impact on the results, albeit we have not observed any significant deterioration of predictions using older versions of the E. coli and B. subtilis genome annotations (data not shown).
We have shown that combining different detection strategies by applying a very simple sequence-dependent filter to our promoter predictions significantly decreases the FP rate. Since accuracy is a trade-off between sensitivity and the FP rate, this procedure could allow the reduction of the threshold, thus leading to the detection of more TP and increasing the sensitivity. The integration of a more sophisticated sequence-dependent method to our strategy could be used to further reduce the FP rate. Distance filters were also successfully used by other groups to decrease the number of FPs [24,39]. However, this can hardly be justified in biological terms as underlined by Huerta and Collado-Vides [24]. Moreover, such filters may not be Bacterial species, gene names, predicted promoter sequences (underlined -35 and -10 boxes), experimentally identified transcriptional start sites (in bold) and the corresponding Z-scores are shown for all tested IRs where a primer extension signal was detected.
suitable for alternative σ factor-dependent promoters or other transcription regulators. We have thus decided not to exploit distance constraints, although it remains possible for an eventual user to determine if a putative promoter is located at an appropriate distance from a gene of interest.

Conclusion
A simple and intuitive concept about the preferential location of regulatory sequences has allowed the identification of principal σ factor dependent-promoter sequences in the genome of various bacteria. Minimal information about the structure of the searched pattern was only required for our algorithm to detect these promoters. Moreover, it could be possible to predict promoters in species for which little transcriptional information is available using the proposed general matrix. Since a biased distribution pattern also appear to be conserved for alternative σ factors and other regulatory proteins in a variety of prokaryotes, it should be possible to design distribution matrices to identify their corresponding DNA binding sites.

Genome-derived distribution matrices
Genomic distributions of hexanucleotide pairs were represented by a ratio of the number of hits in IRs located upstream of a gene (P) to total hits in the whole genome (G). Hits were counted only on the functional strand (on the same strand than the following coding sequence) for all spacings inside the allowed spacing range. Identical hexanucleotide pairs with different spacer length will thus have the same genome-derived distribution matrix provided that their respective spacings are included in the allowed range. The genomic distribution of up to three exclusive mismatches per hexanucleotide was also reported in genome distribution matrices. Ratios at various mismatches combinations were reported in genomederived distribution matrices of dimension 4 × 4 (Figure1, 2). Columns and rows respectively represent mismatches in the -10 (proximal) and -35 (distal) boxes.

Synthetic matrices
248 371 200 distribution matrices were generated in silico and referred to as "synthetic matrices". To create these, the genome-derived distribution matrices of almost all characterized promoter sequences available were analyzed, and the range of variation in each cell was determined in accordance with the observed ratios. The range and step length was independently established in each cell. Detailed information about synthetic matrices is available in Additional file 2.

Score calculation
To calculate a score, the genome-derived distribution matrix of a hexanucleotide pair was compared to a synthetic matrix. The analytical approach was inspired by the image processing field and involved four components, each representing the mean of square differences between matrices: R 1 is calculated on the raw data of the matrices, and R 2 to R 4 are respectively calculated on the horizontal, vertical and diagonal directional derivatives of matrices to evaluate the three different slopes of the matrices. Each slope is related to the representation of the genomic distribution of hexanucleotide pairs (proximal and distal boxes). A weight (w) of ¼ is next applied to each component. The final score = 1/(wR 1 + wR 2 + wR 3 + wR 4 ). See Additional file 3 for a detailed example of score calculation.

IR scanning
IR scanning was accomplished by taking the first six most upstream nt of an enlarged IR along with the downstream hexanucleotide window located at the shortest distance within the specified spacing range. A genome-derived distribution matrix was then generated and a score was calculated with the above described score metric. This procedure was repeated for all allowed spacings by moving the downstream hexanucleotide window by one nt. The upstream hexanucleotide was next moved by one nt and the same procedure was repeated until all appropriate hexanucleotide pairs of the region were processed. The mean of values obtained for all hexanucleotide pairs sharing the same proximal box were then plotted on a graph (Figure 3). Using the maximum values instead of the mean gave very similar results (data not shown). Two thresholds were selected. The region threshold (tR) was set at three standard deviations above the mean of all points from a specific IR. The genome threshold (tG) was set at two standard deviations from the mean of all points from all IR of a genome. tR and tG were optimized using E. coli and B. subtilis promoters data (data not shown). The value of any point had to be higher than both thresholds to be considered as a candidate promoter. All adjacent points above thresholds were combined in one peak and represented by their highest point. The widest peak has 6 points and the mean is 1.25 point per peak. A peak had to be located within 4 nt of an experimentally identified TSS to be considered as a TP. All other points above thresholds were considered as FPs. Points representing characterized hexanucleotide pairs below the highest threshold were considered as FNs, while all other points below this threshold were considered as TNs. According to Tompa et al., sensitivity is defined as TP/(TP + FN), specificity as TN/ (TN + FP), precision (or positive predictive value) as TP/ (TP + FP), and performance as TP/(TP + FN + FP) [38].

Evaluation of synthetic matrices
The evaluation of over 248 million synthetic matrices on the 625 enlarged IRs (containing 684 characterized promoter sequences from the ten genomes mentioned in Table 1) was performed on the Mammouth Linux cluster of the Université de Sherbrooke (1808 processors, 7.6 Tflops, 5.5 TB of RAM, 160 TB of HD). The performance score for each synthetic matrix was calculated and the best matrix selected as the specific matrix for each genome ( Table 1). The general synthetic matrix was selected from the sum of the relative performances of each matrix on each genome following S = (sum i (perf ij /maxPerf j )) where perf ij represents the performance score of a given matrix in the organism j, and maxPerf j represents the maximum performance score of all matrices in the organism j.

Cross-validation tests
Three-fold cross-validation tests were conducted with 1% of the synthetic matrices randomly chosen from the previously described set. The E. coli and B. subtilis datasets were randomly divided into three groups, and all possible combinations of two groups were used to select new specialized synthetic matrices. The statistical indicators were next calculated on the remaining group. Results are presented in Additional file 4.

PSWM scoring function
The scoring function of our method was replaced by a function based on PSWM scores. The rest of the IR scanning procedure remained absolutely identical to the initial design. The promoter datasets of E. coli and B. subtilis were used to construct PSWMs, which were normalized according to the intergenic ATGC content of the tested genome. Results are presented in Additional file 6.

Shuffled genome
Two shuffled genomes were created. First, the regions (genes and IRs) were independently shuffled to conserve the possible AT bias of IR [40]. The second type was made on the entire genome so that no bias is kept. Shuffled genomes were next used to calculate genome-derived distribution matrices to assess the same enlarged IRs previously analyzed (data not shown).

Codon usage evaluation
By definition, a hexanucleotide contains 4 overlapping codons. The mean of the utilization ratio of the eight codons of a hexanucleotide pair was thus compared to the average usage frequency of a codon (15.62/1000 residues for E. coli) to evaluate if there is a difference between hexanucleotide pair sequences precisely identified (TP) or missed (FN) (data not shown). Codon usages were taken from the Codon Usage Database [76].

Biochemical validation of predicted promoters
All IRs for which no promoter sequence is characterized in B. subtilis, E. coli and M. tuberculosis, were analyzed with their respective organism-specialized matrix to predict putative promoters. In order to validate some predicted promoters under the control of the housekeeping σ factor, predictions were selected on the basis of the putative function of their corresponding gene, the Z-score, the loci organization and the promoter sequences. Validation of the M. tuberculosis prediction was made on the closely related non-pathogenic M. bovis BCG-Russia. E. coli K12 ATCC10798 and B. subtilis NIG2001 [77] were grown in LB medium. M. bovis BCG-Russia was grown in Middlebrook 7H9 medium supplemented with Albumine-Dextrose-Saline, Tween 80 and cycloheximide. All cultures were harvested at an OD 600 between 0.6 and 0.8 and RNA was extracted using the Ribopure RNA extraction kit (Ambion) or the RNeasy kit (Qiagen). RNA was quantified by spectrophotometry and integrity was verified on formaldehyde denaturing gel. Primer extensions were performed according to standard procedures. Between 30 and 60 μg of RNA were used for each reaction. Extension products were migrated on 5M urea-6% acrylamide sequencing gels along with sequencing reactions. IRs were cloned in pCR2.1-TOPO TA cloning vector (Invitrogen) or pdrive TA cloning vector (Qiagen). Oligonucleotide primers are listed in Additional file 8. Sequencing ladders were produced with the Sequenase 2.0 kit (USB) according to the manufacturer's instructions. Gels were scanned using a Molecular Dynamics Storm 840 Phosphorimager.


extfieldsofstudy: Biology Medicine

text:
Human Immunodeficiency Virus Persistence and Production in T-Cell Development

ABSTRACT Human immunodeficiency virus type 1 (HIV-1) replication depends on CD4 and coreceptor expression as well as host factors associated with the activation state of the cell. To determine the impact of the activation stage of thymocytes on the HIV-1 life cycle, we investigated R5 and X4 HIV-1 entry, reverse transcription, and expression in discrete thymocyte subsets at different stages of T-cell development. Early after infection, preferential entry and replication of R5 HIV-1 were predominantly detected in mature CD3+/hi CD27+ thymocytes. Thus, R5 HIV-1 targets the stage of development where thymocytes acquire functional responsiveness, which has important implications for HIV pathogenesis. In contrast, X4 HIV-1 expression and replication were primarily found in immature CD3−/+/low CD27− CD69− thymocytes. HIV-1 proviral burden and virus expression in thymocyte subsets correlated with the expression of the highest levels of the respective coreceptor. R5 and X4 HIV-1 entered and completed reverse transcription in all subsets tested, indicating that the activation state of thymocytes and coreceptor expression are sufficient to support full reverse transcription throughout development. Although R5 HIV-1 is expressed mainly in mature CD3+/hi CD27+ thymocytes, 5.3% of HIV-1-infected immature thymocytes express R5 HIV-1, indicating that potentially latent viral DNA can be established early in T-cell development.

Human immunodeficiency virus type 1 (HIV-1) entry involves sequential interactions of the viral envelope proteins with CD4 and chemokine receptors, primarily CCR5 or CXCR4 (11,12,15). With over 90% of developing cells in the thymus expressing CD4, the determining factor for HIV-1 infection in the thymus is the expression of the appropriate chemokine receptor or coreceptor. We and others have described the predominance of CXCR4 expression over CCR5 expression in the thymus, which facilitates the greater entry levels, faster replication kinetics, and enhanced cytopathicity of CXCR4 (X4)-tropic over CCR5 (R5)-tropic viruses in the human thymus (10,18,27,40). Although the levels of CCR5 in the thymus are low, simian immunodeficiency virus studies have shown that the thymus is an important target for R5 virus (17,30,36,49). Sopper et al. found an increase in mature thymocytes in the asymptomatic phase of simian immunodeficiency virus infection, with a severe depletion of mature CD4 and CD8 thymocytes in the symptomatic phase of infection, indicating impairment of T-cell regeneration (36). In addition to entry, the activation state of the host cell is a determining factor for HIV-1 replication.
Cellular activation can either enhance or restrict HIV-1 replication at multiple stages of the viral life cycle (23,45). Completion of reverse transcription in vitro is dependent on the activated state of the host cell at or beyond the G 1 b phase of the cell cycle (21,50,51). Cell-type-specific replication barriers also exist; for example, X4 HIV-1 is reportedly blocked at the level of nuclear import in Th2 T-cell clones (47). Furthermore, a lack of X4 replication and pathogenicity has been observed in the SCID-hu peripheral blood lymphocyte model, in which virus is injected into SCID mice 2 weeks after peripheral blood lymphocyte transfer, when the cells exhibit a memory phenotype (14). Monocytes, macrophages, and dendritic cells also exhibit activation/maturation-related blocks in HIV-1 replication (4,33).
T-cell development in the thymus is an active process of positive and negative selection (37,48). Given the diverse types of signals received by thymocytes during development, specific developmental stages may be more or less permissive for virus replication than others. Thus, HIV-1 pathogenesis in the thymus may be better understood by investigating the replicative capacity of HIV-1 in thymocytes at distinct developmental stages and by examining HIV-1 replication at the activationdependent steps of the viral life cycle within the thymus. Although HIV-1 expression in major thymocyte subsets has been reported, an in-depth analysis of entry and reverse transcription in distinct thymocyte subsets, in particular of thymocytes to be exported to the periphery, is still lacking.
We report here that thymocytes are sufficiently activated to support HIV-1 replication at the five distinct differentiation and activation stages that we examined. Both R5 and X4 HIV-1 entered and completed reverse transcription in each developmental subset studied. However, R5 HIV-1 replication in vivo predominated in the functionally mature CD27 ϩ thymocyte subset while X4-HIV-1 replicated primarily in the CXCR4 ϩ/hi CD69 Ϫ (prepositive-selection) stage of development and was highly productive in the CD3 Ϫ/dim CD71 ϩ cycling cells.
Thymocyte preparation. Normal human postnatal thymus specimens were obtained from children undergoing corrective cardiac surgery. The tissue was placed in NH 4 Cl-Tris lysing buffer to remove the red blood cells while the tissue was cut into small pieces and passed over a cell strainer to generate a single-cell suspension of thymocytes. The cells were washed in phosphate-buffered saline and serum-free medium (AT-IMDM). Thymocytes were resuspended at 2 ϫ 10 7 cells/ml in serum-free medium for HIV-1 infection.
Immunomagnetic selection of thymocyte subsets. Goat anti-mouse M280 magnetic beads (6 ϫ 10 9 to 7 ϫ 10 9 ; Dynal, Lake Success, NY) were diluted to a 60-ml volume with serum-free medium without transferrin (A-IMDM) supplemented with 1% bovine serum albumin. The beads (1 ϫ 10 8 /ml) were coated with 25 l of the following antibodies (50 tests/ml) overnight at 4°C: mouse antihuman CD45RA, mouse anti-human CD69, and mouse anti-human CD71. Thymocyte subsets were selected sequentially as described previously (27). Briefly, total thymocytes were incubated at a 1:1 ratio of CD45RA-coated M280 immunomagnetic beads to cells for 30 to 60 min at 4°C while rotating in a 15-ml culture tube in 4 to 5 ml medium. Cells were separated on a Dynal magnet, and both positively selected bead-bound cells and negatively selected CD45RA Ϫ thymocytes were collected for further rounds of selection. The CD45RA Ϫ population was incubated for 30 min at 4°C with 30 l of CD27 mAb (100 tests/ml) in a 1-ml volume and then washed and incubated with M280 immunomagnetic beads at a 4:1 bead-to-cell ratio for 30 to 60 min at 4°C as described above. Both CD27 ϩ and CD27 Ϫ fractions were collected. The resulting CD45RA Ϫ CD27 Ϫ population was incubated with CD69-coated immunomagnetic beads at a 2:1 ratio as described above. CD69-positive and -negative fractions were collected. The resulting CD45RA Ϫ CD27 Ϫ CD69 Ϫ thymocyte population was incubated with CD71-coated M280 immunomagnetic beads at a bead-to-cell ratio of 1:1 for 30 to 60 min at 4°C while rotating and separated on a Dynal magnet into both positive and negative fractions. Each of the positively selected fractions was subjected to two more rounds of positive selection on the magnet to increase purity. Because the beads could not be removed for flow cytometric analysis of purity, we could only detect the purities of the negative fractions which were between 95 and 99%.
HIV-1 infection of SCID-hu mice and postnatal thymocytes. The syncytiuminducing, CXCR4-tropic hybrid molecular clone HIV-1 NL4-3 (NL4-3) was used for part of these studies (1). Virus stocks were prepared from daily harvests of supernatants from CEM cells (CCRF-CEM) infected with virus derived from COS cells electroporated with plasmid pNL4-3. Stocks of the non-syncytiuminducing, CCR5-tropic molecular clones HIV-1 JR-CSF (JR-CSF) and HIV-1 NFN-SX (NFN-SX) were prepared from 24-h harvests of supernatants from interleukin 2 (IL-2)-and phytohemagglutinin (PHA)-stimulated peripheral blood mononuclear cells infected with the supernatant of COS cells electroporated with plasmids pYKJR-CSF and pSX (22,24). Virus stocks were stored at Ϫ70°C and treated with 2 g/ml DNase (Worthington, Lakewood, N.J.) for 30 min at room temperature in the presence of 0.01 M MgCl 2 before infections. All infections were standardized by quantitating the infectious units in limiting-dilution studies using PHA-stimulated peripheral blood mononuclear cells (50,51).
C.B.17 SCID mice were bred at UCLA and implanted with human fetal thymus and liver (thy/liv) grafts under the murine kidney capsule (2,8). Four to 6 months postimplantation, the thy/liv grafts were infected by direct injection of HIV-1 into the graft. Ten nanograms of p24 from R5 HIV-1 JR-CSF and 2 ng of p24 from X4 HIV-1 NL4-3 were injected into the implants in a 50-l volume. Mock-infected implants, used as controls in all experiments, were prepared by injecting the implants with 50 l of appropriate control supernatants. SCID-hu mice implanted with thy/liv grafts from the same donor were used for each individual experiment. The five different experiments included 6 to 27 SCID-hu mice per experiment, with 6 to 10 mice per time point. Thy/liv graft implants were taken between 1 and 9 weeks postinfection, and thymocytes were isolated as described above for the postnatal thymus specimens.
In vitro infection and culture of thymocytes were performed as previously described (27,43). Briefly, 1 ϫ 10 7 freshly isolated, nonstimulated thymocytes were incubated with 100 ng of X4 and R5 viral p24 in the presence of 10 g/ml Polybrene (Sigma, St. Louis, MO) for 1 h at 37°C. Control thymocytes were mock infected in the presence of Polybrene with supernatants from the same uninfected cells used to prepare the virus stocks. After infection, the cells were washed, resuspended in serum-free medium, and cultured for 18 h.
Intracellular staining for Gag HIV-1 protein was performed after surface immunophenotyping and fixation in 1% paraformaldehyde. Fixed cells were subsequently permeabilized in 0.2% Tween 20 for 15 min at 37°C, washed with phosphate-buffered saline containing 2% newborn calf serum and 0.1% sodium azide (fluorescence-activated cell sorting [FACS] buffer), blocked with human AB serum, and stained with 2.5 l KC57 fluorescent antibody or an IgG control. Finally, cells were washed with 0.2% Tween 20 and resuspended in FACS buffer before acquisition with a dual-laser FACSCalibur flow cytometer (BDIS, San Jose, CA). Multiparameter data acquisition and analysis were performed with Cell Quest software (BDIS).
Cell cycle analysis. Total human thymocytes were surface phenotyped using mAb conjugated to CD71-FITC and CD3-PE, fixed in 0.25% paraformaldehyde for 30 min, permeabilized by 0.2% Tween 20 for 15 min at 37°C, and incubated with 20 ng/ml of the DNA dye 7-amino-actinomycin D. Samples were acquired with a FACScan flow cytometer (BDIS, San Jose, CA). The FL3 area and FL3 width plot were displayed to gate out any aggregates with more then 4N DNA and subapoptotic cells with less then 2N DNA content. Data were analyzed using the cell cycle analysis software ModFit LT (from Verity Software House, Inc., Topsham, ME).
DNA isolation and PCR. DNA was isolated from in vitro-infected cultures at 18 h postinfection by phenol-chloroform extraction and ethanol precipitation as described previously (27). For samples collected from SCID-hu thymocyte subsets, simultaneous extraction of RNA and DNA was performed using TRIzol (Gibco/BRL-Life Technologies), following the manufacturer's protocol. Quantitative PCR was performed using an ABI Prism 7700 sequence detection system in a 25-l reaction volume for 40 cycles between 95°C for 15 seconds and 60°C for 1 min, using the components in the TaqMan core reagent kit (Applied Biosystems, Foster City, CA). Proviral burden was measured using the primer pair SR1/AA55, which binds within the R and U5 regions of the long terminal repeat (LTR). This same pair was used to detect nucleotides synthesized at the initiation of reverse transcription. The SR1/661 primer pair spans the primer binding site and measures the completion of reverse transcription. The primers SR1 (5Ј-CAAGTAGTGTGTGCCCGTCTGT-3Ј), AA55 (5Ј-CTGCTAGAGAT TTTCCACACTGAC-3Ј), and 661 (5Ј-CCTGCGTCGAGAGAGCTCCTCTG G-3Ј) were used along with the fluorogenic probe ZXF (5Ј-TGTGACTCTGGT AACTAGAGATCCCTCAGACCC-3Ј), modified with 6-carboxyfluorescein reporter dye on the 5Ј end and 6-carboxytetramethylrhodamine quencher dye on the 3Ј end. To control for the amount of cellular DNA per reaction, the following primers and probe were used to amplify a portion of the human beta-globin gene: forward primer BGF1 (5Ј-CAACCTCAAACAGACACCATGG-3Ј), reverse primer BGR1 (5Ј-TCCACGTTCACCTTGCCC-3Ј), and the fluorogenic probe modified as described above (5Ј-CTCCTGAGGAGAAGTCTGCCGTTACTG CC-3Ј).
Statistics. The Wilcoxon signed rank test was used to compare differences in HIV-1 entry, proviral burden, and levels of coreceptor expression for thymocytes at different stages of maturation. P values were adjusted for intracorrelation and intraexperiment correlations, where appropriate, through permutation methods. Due to small sample sizes, attention has been drawn to results that were significant at the 0.10 level or lower.

RESULTS
Delineation of distinct developmental stages of thymopoiesis. As T cells develop in the thymus, they receive developmental cues driving T-cell receptor (TCR) gene rearrangement, expression, and selection. These developmental stages are identified by the progressive expression of cell surface markers, leading to the eventual phenotype of naïve CD4 and CD8 single-positive T cells in the periphery (Fig. 1A). Using four-color flow cytometry, we immunophenotyped thymocytes with various combinations of the following activation/developmental stage markers: CD1, CD3, CD4, CD8, CD45RA, CD27, CD69, and CD71 (for our purposes indicated as stages I to IV and substage I-P). The most mature thymocytes at the last stage of thymopoiesis, about to emigrate from the thymus, express high levels of CD3 and the leukocyte common antigen CD45RA isoform (stage IV) (Fig. 1A). These mature CD45RA ϩ CD3 ϩ/hi thymocytes express either CD4 or CD8 (16,41,44) (Fig. 1A). Thymocytes at the developmental stage preceding CD45RA expression (CD45RA Ϫ CD3 ϩ/hi ) express the tumor necrosis factor receptor family member CD27 (stage III) and include late CD4/CD8 double-positive cells as well as CD4 and CD8 single-positive cells (Fig. 1A). CD27 has been shown to function as a costimulatory molecule affecting the proliferation of peripheral T cells, which is consistent with CD27 expression in the thymus, where CD27 expression marks thymocytes that gain functional proliferative responsiveness to mitogens such as PHA (28,46). The more immature population preceding CD3 ϩ/low CD45RA Ϫ CD27 Ϫ cells (stage II) bears the activation marker CD69, a C-type lectin binding protein that arises in response to positive selection signaling for an appropriate TCR rearrangement (42). The most immature cells we examined were cells at the stage before positiveselection CD3 Ϫ CD45RA Ϫ CD27 Ϫ CD69 Ϫ cells (stage I).
Within this most immature stage, a subset of cells expresses the transferrin receptor CD71, also referred to as the proliferating antigen (stage I-P). In the thymus, the CD71 ϩ population is enriched in cells with more than 2N DNA as determined by flow cytometric measurements of DNA content by the DNA fluorescent molecule 7-amino-actinomycin D (Fig. 1B).
X4 and R5 HIV-1 replicate at distinct stages of thymopoiesis in vivo. To further delineate the stages of T-cell activation or development at which R5 and X4 HIV-1 replicate, thy/liv graft implants in SCID-hu mice were infected with X4 HIV-1 (NL4-3) (four experiments) or R5 HIV-1 (JR-CSF) (four experiments). Despite the rapid replication kinetics of NL4-3, virus expression was not detectable at 1 week postinfection in five out of five thy/liv grafts studied at that time point. At 19 days after infection with NL4-3, virus expression was detected by intracellular staining for Gag HIV-1 in two out of two thy/liv grafts at frequencies of 0.05 to 0.3% of the total thymocytes. At 4 and 5 weeks after infection with NL4-3, 6 out of 10 of the thy/liv graft implants expressed virus in 5 to 9% of the total thymocytes and depletion of CD4 ϩ and CD4 ϩ CD8 ϩ cells was observed. To detect JR-CSF expression, infected thy/liv graft implants were collected at weeks 1, 3.5, 4, 7, 8, and 9. No virus expression was detected at 1 week postinfection in five out of five JR-CSF-infected mice. However, 31 out of 33 thy/liv graft implants expressed JR-CSF after week 3.5 postinfection and the numbers of Gag HIV-1 ϩ thymocytes ranged from 0.15 to 1.5% of the total population. By intracellular staining for Gag HIV-1 protein, we calculated the frequency of virus-expressing cells at each of the developmental stages of thymocyte maturation. In NL4-3-infected implants, more than 80% of the virus-expressing cells were present in the immature CD69 Ϫ (stage I; CD3 Ϫ/low ) popula- Ϫ thymocytes; stage I-P (a subset of stage I in which thymocytes express CD71), CD69 Ϫ thymocytes; stage II, CD69 ϩ cells that have not yet acquired CD27 expression; stage III, CD69 ϩ CD27 ϩ thymocytes that have not yet acquired CD45RA expression; stage IV, CD27 ϩ thymocytes expressing CD45RA. Within stage I, the first circle represents the CD34 ϩ thymic immigrant. The second circle represents the immature CD4 single-positive thymocyte. The third and fourth circles represent CD3 Ϫ CD4 ϩ CD8 ϩ doublepositive (DP) thymocytes and CD3 ϩ/low DP thymocytes, respectively. The fifth circle (stages II/III) represents CD69 ϩ positively selected thymocytes in transition to the next stage of either CD4 single-positive or CD8 single-positive thymocytes. The last stage (stage IV) represents the CD3 ϩ/hi CD45RA ϩ thymic emigrant. (B) Cell cycle analysis was determined for total and CD3 Ϫ CD71 ϩ thymocytes by a combination of cell surface staining with antibodies to CD3 and CD71 and 7-AAD. Histograms were drawn for DNA content (7-AAD) in ungated thymocytes (left panel) and CD3 Ϫ CD71 ϩ gated thymocytes (right panel). The proportions of the thymocytes in the G 0 /G 1 stages and S/G 2 /M stages of the cell cycle were calculated by Modfit software.
VOL. 13,2006 HIV-1 REPLICATION IN DISTINCT T-CELL DEVELOPMENT STAGES 1239 tion (Fig. 2). In contrast, in JR-CSF-infected SCID-hu mice, more than 70% of the virus-expressing thymocytes were mature CD3 ϩ/hi cells (Fig. 2). The majority of these mature JR-CSF-expressing thymocytes were CD27 ϩ and CD69 ϩ , and 24% expressed CD45RA. Thus, at early time points after in vivo HIV-1 infection, X4 HIV-1 was expressed in very immature (stage I) thymocytes, while R5 HIV-1 was expressed mainly in mature (stage III) thymocytes. Differences in R5 and X4 HIV-1 proviral burden in thymocyte subsets. To determine whether the differential expression of HIV-1 by X4 and R5 HIV-1 in thymocyte subsets was due to differences in the proviral levels or to a block in virus expression, we examined the HIV-1 DNA content in the thymocyte subsets by quantitative DNA PCR. At 4 and 7 weeks after infection of the thymic implant, we analyzed the R5 HIV-1established infection in each of the thymocyte subsets. The proviral levels were calculated as numbers of copies per 100 cell equivalents as follows: for stage I-P, 10 Ϯ 3.7; for stage I, 8.7 Ϯ 2.8; for stage II, 5.1 Ϯ 1.6; for stage III, 13.4 Ϯ 3.8; for stage IV, 8.6 Ϯ 0.7; and for the total population, 6.9 Ϯ 2.5 (Fig.  3). These data indicate that proviral DNA was also high in immature thymocyte subsets. However, the R5 proviral burden in the stage III thymocyte subset was significantly higher (P ϭ 0.04) than those in other subsets.
The average percentage of thymocytes expressing R5 Gag HIV-1 /cell (Table 1) was less than the average number of proviruses detected in this population (assuming one virus infection per cell in each subset) (Fig. 3), indicating that R5 HIV-1 persists through development in a nonexpressed form and is present as early as stage I-P (in the CD3 Ϫ /CD71 ϩ population). Dividing the percentage of cells expressing viral proteins (Table 1) by the percentage of provirus-harboring cells (Fig. 3), we found that the percentages of cells in each subset expressing R5 HIV-1 are as follows: for stage I-P, 12%; for stage I, 5.3%; for stage II, 17.8%; for stage III, 19.2%; and for stage IV, 26.7%.
As X4 HIV-1 causes substantial CD4 ϩ thymocyte depletion at weeks 4 to 7 postinfection, we chose to analyze X4 HIV-1 infection before T-cell loss occurs. We therefore determined the X4 proviral burden early, at 19 days postinfection. The proviral burdens, calculated as numbers of copies per 100 cell equivalents, were as follows: for subset I-P, 21.3 Ϯ 9.8; for subset I, 8.5 Ϯ 6.7; for subset II, 3.9 Ϯ 3.5; for subset III, 2.3 Ϯ 2.1; for subset IV, 1.1 Ϯ 1.2; and for the total population, 0.5 Ϯ 0.6. The X4 proviral burden in the stage I-P thymocyte subset was significantly higher (P ϭ 0.07) than those in other subsets. We found that the subset producing most of the X4 virus also harbored the greatest proviral burden (Table 1 and Fig. 3).
Presumably, nonexpressed X4 HIV-1 is established at early developmental stages as well; however, with such rapid replication and CD4 thymocyte destruction in the X4 HIV-1 NL4-3infected implant, not enough cells were available to select and purify each subset at late stages postinfection. Certainly, by 19 days postinfection, high X4 HIV-1 burden was achieved (Fig.  3), with barely detectable expression in the CD3 Ϫ CD71 ϩ thymocyte subset (Table 1). Thus, all through thymocyte de- velopment, cellular conditions may lead to a certain degree of nonexpressed R5 and X4 viral DNA or potentially latent HIV-1. CCR5 and CXCR4 expression in distinct thymocyte subsets. Four-color flow cytometric analysis of CXCR4 and CCR5 expression was performed among the five stages of activation presented above to compare the relative levels of coreceptor expression within the thymocyte subsets. As we and others have described previously, the percentages of cells expressing CCR5 were low at all stages, with mature medullary thymocytes expressing the highest percentage of CCR5 (6,26). We show more specifically that the greatest expression of CCR5 was present at stage III (CD27 ϩ ), with 2.5% Ϯ 0.8% (Fig. 4). The percentage of CCR5 expression for stage III thymocytes was significantly higher (P ϭ 0.04) than those for thymocytes at other stages of differentiation. In contrast, CXCR4 was expressed at high percentages in all thymocyte subsets. The highest levels of CXCR4 expression, as measured by mean fluorescent intensity, occurred on immature thymocytes and particularly on immature thymocytes at the cycling stage (stage I-P; CD45RA Ϫ CD27 Ϫ CD69 Ϫ CD71 ϩ ) (Fig. 4) (6,20,27). The level of CXCR4 expression on stage I-P thymocytes was significantly higher (P ϭ 0.07) than those on thymocytes at other stages of differentiation. At the point of positive selection (stage II), CXCR4 was downregulated (39) (Fig. 4). Finally, prior to exiting the thymus, thymocytes demonstrated increased expression of CXCR4, and the frequency of CCR5 expression declined, generating the coreceptor levels found for peripheral blood naïve T cells (Fig. 4). Thus, differential levels of HIV-1 coreceptor expression were observed in thymocyte subsets at different stages of maturation. These patterns of coreceptor expression were found in all five postnatal thymus tissues examined.
HIV-1 entry in discrete thymocyte subsets correlates with levels of coreceptor expression. HIV-1 provirus was detected in all thymocyte subsets tested, yet the virus may have entered at an earlier developmental stage. Therefore, to examine the susceptibility of each of the thymocyte subsets to direct HIV-1 infection, we purified each of the five subsets (stages I-P and I to IV) and infected each subset in vitro. The role of Env for conferring the chemokine receptor tropism of a virus was examined by infection of thymocyte subsets with R5 HIV-1 NFN-SX , derived from NL4-3, with only 157 amino acid differences in the Env protein (24). Each of the thymocyte subsets was infected in parallel with X4 HIV-1 NL4-3 , R5 HIV-1 NFN-SX , and R5 HIV-1 JR-CSF in vitro for 18 h. As we and others described previously, X4 entry was greater than R5 entry in the thymocyte subsets due to the greater prevalence of CXCR4 expression than CCR5 expression in thymocytes (6,27). Viral entry of R5 HIV-1 (NFN-SX and JR-CSF) was highest in the CD27 ϩ CD45RA Ϫ thymocyte subset in three out of three experiments in which cells were infected with equal amounts of virus (Fig.  5). The mean viral-copy numbers of NFN-SX and JR-CSF a Thymocytes were isolated from HIV-1-infected implants from 2.5 to 7 weeks postinfection. The cells were surface stained with antibodies to CD3, CD71, CD69, CD27, and CD45RA combined with intracellular staining for the Gag HIV-1 protein (KC57-FITC). Each specific subset was gated electronically, and the frequencies of Gag HIV-1 (KC57 ϩ )-expressing cells in the subsets were determined. ND, not determined. DNA per 1,000 cells were as follows, respectively: for the total thymocytes, 2 and 2.4; for stage I-P, 10.5 and 1.5; for stage I, 3.3 and 3.9; for stage II, 4.2 and 3; for stage III (CD27 ϩ ), 23 and 7.6; and for stage IV, 6 and 5.2. These data indicate that R5 HIV-1 viral entry was 2-to 12-fold higher in CD27 ϩ cells than in the total population after in vitro infection ( Fig. 5A and B). R5 HIV-1 NFN-SX entry was statistically significantly higher in stage III thymocytes (P ϭ 0.02) than in other subsets, while HIV-1 JR-CSF entry in stage III thymocytes reached statistical significance (P ϭ 0.09). R5 HIV-1 JR-CSF entry was dependent on the presence of CCR5 since entry was prevented by the CCR5 antagonist TAK779 (3) (data not shown). X4 HIV-1 entry was quite efficient among all thymocyte subsets but was significantly higher in stage I-P thymocytes (P ϭ 0.01) than in other subsets. Thus, X4 and R5 HIV-1 entered thymocyte subsets at all stages of differentiation. However, while R5 entry occurred in all subsets, the highest level of R5 entry was observed in the CD27 ϩ CD45RA Ϫ developmental stage (stage III). Each thymocyte subset is sufficiently active for R5 and X4 HIV-1 reverse transcription. As reported for peripheral quiescent T cells, the completion of reverse transcription is dependent on the activation state of the cell (21). By harvesting FIG. 4. HIV-1 coreceptor expression in thymocytes at five developmental stages. CCR5 and CXCR4 expression at stages I-P to IV of thymocyte development was determined by four-color flow cytometric analysis using combinations of monoclonal antibodies to CCR5-PE, CXCR4-PE, CD3-FITC and -APC, CD45RA-APC, CD27-FITC, CD69-TC, and CD71-FITC in five separate experiments as described in Materials and Methods. Asterisks indicate statistically significant differences. As only small numbers of thymocytes express CCR5 and almost all thymocytes express CXCR4, CCR5 data are expressed in percentages and CXCR4 data as mean fluorescence intensities. (A) The mean percent CCR5 ϩ cells in each subset was as follows: for subset I-P, 0.95% Ϯ 0.6%; for subset I, 0.08% Ϯ 0.07%; for subset II, 0.16% Ϯ 0.06%; for subset III, 2.5% Ϯ 0.8%; and for subset IV, 1% Ϯ 0.9%. (B) Since CXCR4 is expressed in each of these subsets, we measured the mean fluorescent intensity (or geometric mean) of CXCR4 for each subset relative to that for subset I-P, which consistently expressed the highest levels of CXCR4 in all five experiments. Subset I-P, 100%; subset I, 32% Ϯ 9%; subset II, 14% Ϯ 9%; subset III, 23% Ϯ 17%; and subset IV, 39% Ϯ 28%.

DISCUSSION
Upon examining HIV-1 replication in thymocyte subsets associated with markers of activation, we found for the first time preferential entry and replication of R5 HIV-1 more specifically in CD3 ϩ/hi CD27 ϩ thymocytes in vivo in the SCID-hu model. This finding has important implications for HIV-1 pathogenesis, as thymocytes acquire functional capabilities at the CD3 ϩ/hi CD27 ϩ CD1 Ϫ stage of development (28). Thus, despite the relatively low levels of viral replication and depletion of CD4 ϩ thymocytes after R5 HIV-1 infection, infection of a functionally important subset can potentially severely impact the regeneration of functional T cells in the periphery. We further show that the highest level of X4 HIV-1 expression occurred in the immature CXCR4 ϩ/hi CD3 Ϫ CD71 ϩ proliferating thymocytes at early time points after infection. These data are in agreement with previous reports that R5 HIV-1 replicates in the mature thymocytes in the medulla (5,20,29) and X4 HIV-1 replicates in the immature, CXCR4 ϩ/hi thymocytes in the cortex of the thymus but further delineates the precise stage of thymopoiesis for R5 and X4 HIV-1 targets.
The CD3 ϩ/hi CD27 ϩ CD1 Ϫ thymocyte subset marks cells that have received signals through the T-cell receptor to drive positive selection of the CD4 ϩ CD8 ϩ thymocytes and have acquired functional responsiveness (28). Thus, after positive selection, CXCR4 expression decreases and CCR5 expression increases. This is also the stage where thymocytes respond to the cytokines IL-2 and IL-4 to produce high levels of HIV-1 in vitro (26). This process of decreased CXCR4 and increased CCR5 expression after T-cell receptor signaling in the thymus mirrors an occurrence in the periphery, where CXCR4 expression decreases and CCR5 expression increases as naïve cells progress to memory cells (7). In addition, functional responsiveness to antigen is improved in the memory cells compared to that in the naïve T cells, just as functional responsiveness to mitogens is acquired in the positively selected thymocytes (28).
Despite the low level of CCR5 expression, we found that R5 viruses could infect and replicate in all the subsets tested. In agreement with our results, Scoggins et al. showed that even a more cytopathic R5 HIV-1 isolate persisted in using exclusively CCR5 in the SCID-hu model (35). As predicted by CXCR4 expression, X4 viruses also entered all the subsets tested; however, the greatest viral burden and highest frequency of expression were observed in the CD3 Ϫ CD71 ϩ subset. Consistent with CD71, the transferrin receptor, which is referred to as the proliferation antigen, we found about 50% of the CD71 ϩ cells well into the S/G 2 and M phases of the cell cycle, whereas only 20% of the total population was beyond the G 1 phase of the cell cycle. Thus, part of the enhanced cytopathicity of X4 viruses in the thymus may be due not only to direct cell depletion but also to the destruction of the replicative capacity of early thymocytes. Camerini et al. calculated a proviral load of one provirus for 11 cells as a threshold for CD4 ϩ thymocyte depletion in the SCID-hu mouse model (10). By our measurements, X4 HIV-1 NL4-3 achieved this threshold in immature thymocytes at the CD3 Ϫ CD71 ϩ stage (stage I-P) as early as 19 days postinfection, while R5 HIV-1 reached this threshold at the mature CD3 ϩ/hi CD27 ϩ stage (stage III).
We also show that human thymocytes at multiple stages of development are permissible for direct HIV-1 entry and reverse transcription in vitro. Though quiescent peripheral T lymphocytes render HIV-1 unable to complete reverse transcription in vitro (50), we found that thymocyte subsets at numerous developmental stages are sufficiently activated for HIV-1 to complete reverse transcription in vitro. Factors such as TCR signal strength and IL-2 signaling of resting peripheral T cells are needed for HIV-1 infection to pass the reverse transcription stage (25). Thymocytes may therefore receive sufficient signals by endogenous cytokines or TCR engagement for completion of reverse transcription to take place. Alternatively, Korin and Zack indicated that completion of reverse transcription was dependent on peripheral T cells at or beyond the G 1 b phase of the cell cycle (21). Though we did not determine the cell cycle stages of the various thymocyte subsets, we did not observe any limitations in the reverse transcription process for any of the subsets tested. However, a block in vitro may not correlate with a block in vivo, for Eckstein et al.
showed that in vitro nonproductively infected naïve quiescent T cells are productively infected when residing in the microenvironment of the lymph nodes in vivo (13). Thus, HIV-1 replication in distinct thymocyte subsets correlates more with coreceptor expression than with the activation state of thymocytes.
Human thymocytes appear permissible to HIV-1 replication by well-characterized molecular clones NL4-3 and JR-CSF. However, Stoddart et al. reported that protease inhibitor-resistant isolates are not propagated in the SCID-hu thymus despite sufficient replication in peripheral T cells (38). The mechanism for this thymocyte-selective block in HIV-1 replication is unknown, but it does raise the question as to how well molecularly disparate quasi-species present in vivo are able to replicate in the thymus. We previously measured and compared the replicative capacities of 10 HIV-1 isolates from children to those of the well-characterized X4 HIV-1 NL4-3 and R5 HIV-1 JR-CSF clones in vitro (26). R5 HIV-1 isolates prevailed over X4 isolates among these pediatric isolates, which were obtained as close to birth as possible. We found that X4 primary isolates were not more cytopathic for thymocytes than the R5 isolates, indicating that coreceptor use by HIV-1 primary isolates was important but not the sole determinant of HIV-1 pathogenesis in the thymus (26). These data are supported by the finding that replication in CD4 ϩ T-cell cultures of R5 HIV-1 isolates is higher than that in X4 and dual tropic HIV-1 isolates (34). Camerini et al. reported that the determinants for cytopathicity were due to replication levels, and X4 Env enabled greater entry and replication than R5 HIV-1 in the thymus, reaching the critical threshold of replication for cytopathic effects (10).
Brooks et al. proposed that X4 HIV-1 latency in peripheral T cells could be established after infection of the thymus by using the SCID-hu mouse model (9). Latency formation is attributed to a decrease in viral transcription as thymocytes mature. In agreement with this, we found that the frequencies of X4-expressing cells decreased within progressively maturing thymocyte subsets. However, with R5 HIV-1, we found that the frequencies of R5 virus expression and proviral burden increased in thymocytes at more-mature stages of development (stages III and IV). Although R5 HIV-1 replication frequency increased in the mature subsets, at stages III and IV, where 2 to 5% of the cells expressed HIV-1, the proviral burdens ranged from 8 to 20% (on average, greater than 80% of R5 HIV-1-infected thymocytes did not express HIV-1), thereby confirming the notion that a pool of potentially latently infected cells is likely exported to the periphery. Even at the early (I and I-P) stages of thymopoiesis, R5 HIV-1 expression occurred in 0.3 to 3% of thymocytes, whereas the proviral burdens ranged from 10 to 28% of thymocytes in stages I and I-P, indicating that latency can be potentially established early in development and maintained throughout thymocyte development. Thus, thymocytes at multiple stages of development are permissive for viral replication and are at an optimal activation state for viral entry, reverse transcription, and latency formation.
Our results were obtained with the SCID-hu mouse model, which allows access to, and dissection of, thymocyte subsets for extensive ex vivo analysis of R5 and X4 HIV-1 infection in thymocyte subsets. The model was the first to correctly show that nef deletion results in attenuated HIV-1 replication, that breakthrough of HIV-1 replication can occur in vivo in the face of antiretroviral therapy, and that the thymus is not permanently damaged by HIV-1 replication (reviewed in reference 19). These findings have consistently shown that the SCID-hu mouse model accurately recapitulates HIV-1 pathogenesis in humans. There are certainly limitations to the model, such as a lack of immune responses and an absence of peripheral lymphoid cells (19). However, the SCID-hu mouse model does allow for an accurate assessment of the impact of HIV-1 on the development and regeneration of naïve T cells.


extfieldsofstudy: Medicine Physics Biology

text:
An Ising-Like model for protein mechanical unfolding

The mechanical unfolding of proteins is investigated by extending the Wako-Saito-Munoz-Eaton model, a simplified protein model with binary degrees of freedom, which has proved successful in describing the kinetics of protein folding. Such a model is generalized by including the effect of an external force, and its thermodynamics turns out to be exactly solvable. We consider two molecules, the 27th immunoglobulin domain of titin and protein PIN1. In the case of titin we determine equilibrium force-extension curves and study nonequilibrium phenomena in the frameworks of dynamic loading and force clamp protocols, verifying theoretical laws and finding the position of the kinetic barrier which hinders the unfolding of the molecule. The PIN1 molecule is used to check the possibility of computing the free energy landscape as a function of the molecule length by means of an extended form of the Jarzynski equality.

The mechanical unfolding of proteins is investigated by extending the Wako-Saitô-Muñoz-Eaton model, a simplified protein model with binary degrees of freedom, which has proved successful in describing the kinetics of protein folding. Such a model is generalized by including the effect of an external force, and its thermodynamics turns out to be exactly solvable. We consider two molecules, the 27th immunoglobulin domain of titin and protein PIN1. In the case of titin we determine equilibrium force-extension curves and study nonequilibrium phenomena in the frameworks of dynamic loading and force clamp protocols, verifying theoretical laws and finding the position of the kinetic barrier which hinders the unfolding of the molecule. The PIN1 molecule is used to check the possibility of computing the free energy landscape as a function of the molecule length by means of an extended form of the Jarzynski equality. Manipulation experiments on single biomolecules have greatly increased our knowledge of the structural properties of such molecules. In a typical experiment a controlled force is applied to one of the free ends of the molecule, and the induced elongation is measured. Such experimental techniques have been used to probe the structure of proteins [1,2] and nucleic acids [3,4]. According to the common interpretation the unfolding of a molecule being pulled from one of its ends is hindered by kinetic barriers associated with the strongest linkages which serve to stabilize the molecular structure. The breaking of a molecular bond can thus be viewed as the overcoming of a kinetic barrier. It has been argued [5] that the study of the kinetics of bond breaking under different loading rates can provide much information about the internal structure of molecules, and in particular allows one to measure the strength of the molecular bonds, and to associate to them a position along the molecular structure. In the case of simple molecules, such as RNA hairpins, it is easy to obtain information on the molecular structure by pulling experiments [3]. However, when one deals with large molecular structures, such as multi-domain proteins, the inference of structural characteristics from the unfolding kinetics can prove a difficult task. Therefore, the study of simple models for the unfolding of proteins, whose microscopic native structures are known a priori, is highly desirable: investigating the kinetics of such models can shed new light on the relations between the experimentally observed unfolding features and the molecular structures.
The experiments discussed above are usually performed in non-equilibrium conditions: because of technical limitations the pulling process is faster than the typical molecular relaxation time. The problem of irreversibility of unfolding processes can be avoided by using the remarkable equality introduced by Jarzynski [6], which allows one to measure the free energy difference between the folded and the unfolded state of a biomolecule [7]. By using an extended form of the Jarzynski equality (JE) the free energy landscape of simple models of biopolymers has been probed as a function of the molecular elongation [8,9,10]. Although this approach appears very promising, it still has to be tested on systems of increasing complexity.
Here we approach the mechanical unfolding problem by means of a suitable generalization of the Wako-Saitô-Muñoz-Eaton (WSME) protein folding model [11,12,13]. This is a simplified statistical mechanical model where a binary variable is associated to each peptide bond. The equilibrium thermodynamics has been solved exactly [14,15] and the model has been quite successful in describing the kinetics of protein folding [16,17,18,19] and has also found applications in different fields (see [20] and references therein).
In the present paper we first extend the WSME model by considering the effect of an external force, and show that the equilibrium properties of this new model can be computed exactly, similarly to the original WSME model. In order to mimic the mechanical unfolding of proteins, we use computer simulations and study the unfolding kinetics of our model, both in the cases of constant force and dynamic loading. Finally we probe the free energy landscape of a model protein, exploiting the JE.
The WSME model describes a protein of N + 1 aminoacids as a chain of N peptide bonds (connecting consecutive aminoacids) that can live in two states (native and unfolded) and can interact only if they are in contact in the native structure and if all bonds in the chain between them are native. To each bond is associated a binary variable m k , with values 0, 1 for unfolded and native state respectively. The effective Hamiltonian of the model reads (1) where k B is the Boltzmann constant, and T is the absolute temperature. The first term assigns an energy ǫ i j < 0 to the contact (defined as in [13,14]) between bonds i and j if this takes place in the native structure (∆ i j = 1 in this case and ∆ i j = 0 otherwise). The second term represents the entropic cost q i > 0 of ordering bond i. In Ref. [14] it is shown how to compute exactly the partition function Z = {m k } exp[−βH 0 ({m k })] and the corresponding thermal averages. Here and in the following, the quantity β indicates the inverse temperature β = 1/(k B T ).
In order to couple the protein to an external force we assume that a configuration {m k } of the model defines a sequence of native stretches, separated by unfolded peptide bonds (see inset of Fig. 1). Peptide bonds i and j delimit a native stretch if and only if m i = m j = 0 and m k = 1 for i < k < j. A native stretch is regarded as a rigid portion of the molecule, even under application of the external force, with an end-toend length l i j . If j = i + 1 the stretch reduces to a single aminoacid, which is also regarded as a rigid structure, with length l i,i+i . The values of the parameters l i j are taken from the native structure of the protein [21]. Boundary conditions are introduced through the dummy bonds m 0 = m N+1 = 0. Given the direction of the external force, we assume that a native stretch, or a single aminoacid delimited by two successive unfolded bonds, can only take two orientations, parallel or antiparallel to the force, so that they contribute ±l i j to the length of the molecule. Therefore, given a configuration {m k } of the peptide bonds, we introduce a variable σ i j for each rigid portion of the molecule, either a native stretch or an aminoacid delimited by two successive unfolded bonds, taking values +1 (rigid portion parallel to the force) or −1 (antiparallel). Thus the end-to-end length of the molecule, in the force direction, reads Let us define the Hamiltonian H as the sum of the interaction energy term, contained in the effective Hamiltonian H 0 (1), and of the term − f L, which takes into account the effect of the external force Note that the definition of molecular length (2) is such that the set of variables {σ i j } is dynamically defined by the bond configuration {m k }: for each configuration {m k }, we consider only those variables σ i j such that and, as a function of {m k }, has the same structure as the WSME model (see Eq. (1)). Therefore its equilibrium thermodynamics is exactly solvable, as discussed in Ref. [14]. In the case f = 0, H eff reduces to Eq. (1) (up to an additive constant), with q i = ln 2. Thus, the effective entropic terms can be viewed as resulting from microscopic orientational degrees of freedom of the native stretches.
In the following, the quantity ǫ will indicate the system energy scale. This quantity, together with the interaction energies ǫ i j , will be chosen as in [13,14,18], see also [21]. We also introduce the reduced temperature in terms of the energy scale: T = k B T/ǫ . The time scale will be indicated by t 0 .
Following Ref. [14], for any choice of the model parameters, one can calculate the equilibrium value of the order parameter m = 1/N k m k and of the molecule length, as defined by Eq. (2), for varying force and temperature. We first consider the titin immunoglobulin domain I27 (89 aminoacids, pdb code 1TIT), which has been widely studied both experimentally [1] and theoretically [22]. The energy scale for such a molecule is taken to be ǫ/k B = 43.1 K, while the melting temperature is T m ≃ 8.03 [21]. In figure 1, the root mean square length of the 1TIT molecule is plotted as a function of the external force f for different temperatures. The plateau appearing in Fig. 1, in the low temperature regime, corresponds to the overall alignment of the molecule in the native configuration (m ≃ 1) to the force direction. Having introduced the molecular length we have built a framework within which the mechanically induced protein unfolding can be simulated by applying an external force. In the following, two manipulation schemes will be considered: in the first one the molecule is manipulated in a "force-clamp", where a sudden force is applied to one of the molecule's free ends, in the second one a time-varying force is applied, so that the load on the molecule increases gradually. In order to study the molecular unfolding, we run Monte Carlo simulations with Metropolis kinetics using the Hamiltonian H({m k } , {σ i j }, f (t)). In the following the time scale t 0 will correspond to a single Monte Carlo step. In the force clamp manipulation experiments, the molecule unfolds after a given time τ u which fluctuates between one realization of the unfolding process and the other, due to the stochasticity of the unfolding process. Unfolding can be viewed as an activated process [23], whose kinetics is dominated by a characteristic energy barrier ∆E u placed at the value x u of the reaction coordinate: therefore, it is usu-ally assumed that the average unfolding time τ u follows an Arrhenius law τ u = ω −1 0 exp β (∆E u − f x u ) , where ω 0 is a characteristic attempt rate depending on the microscopic features of the system.
We simulate force clamp manipulations of the 1TIT molecule as follows. Starting from thermal equilibrium with f = 0, at t = 0 we apply a non-zero force f and thus measure the unfolding time τ u as the first passage time of the molecule length across the threshold value L u , defined as half the molecule maximal length, L u = 140 Å, see figure 1. In figure 2 the unfolding time τ u of the 1TIT molecule, averaged over 1000 independent unfolding trajectories, is plotted as a function of the applied force f , for three values of the temperature T = 4, 6, 8. From a fit of the data shown in Fig. (2) to the Arrhenius law, we find that the unfolding length is x u ≃ 3 Å, and is independent of the temperature, as expected (see caption of fig. (2) for the exact values and uncertainties). Note that this value is in good agreement with the experimentally measured value of the titin unfolding length x u = 2.5 Å, found in Ref. [1]. We now consider the case where a time-dependent force is applied to our model molecule and the unfolding time is sampled over 1000 independent trajectories. Here the force increases linearly with time, with a rate r, and thus the rupture force f * is given by f * = rτ u , where unfolding time is defined as in the case of the force clamp. In refs. [23] it has been argued that, if the energy barrier ∆E u is large (compared to the thermal energy k B T ) and rebinding is negligible, the typical unbinding force of a single molecular bond under dynamic loading is given by where τ 0 is the characteristic unfolding time at zero force, τ 0 = ω −1 0 exp β∆E u . In Fig. 3 the breaking force f * is plotted as a function of the pulling velocity, for the 1TIT molecule, for three values of the temperature. The value of the unfolding length obtained by fitting the data to Eq. (4) is x u ≃ 3 Å, and is independent of the temperature, as expected (see caption of fig. 3 for the exact values). Note that this value agrees with that found with the force clamp manipulation, and with the experimental value x u = 2.5 Å found in Ref. [1]. On the other hand, in recent works [24], it has been argued that the rupture force f * has a more complex where the exponent ν depends on the microscopic details of the energy landscape, and γ is the Euler-Mascheroni constant γ ≃ 0.577. This equation reduces to Eq. (4) in the case ν = 1 or in the limit ∆E u → ∞ [24]. In the inset of Fig. 3, we plot the fits of the rupture force data to the equation definingf * . Although the agreement of the data with the equation appears to be rather good, the statistical errors of the fit parameters are quite large, sincef * depends nonlinearly on the set of the unknown parameters. This issue will be addressed in a forthcoming paper. We now evaluate the effective free energy landscape as a function of the molecular length L of the extended model here discussed. Formally, the free energy function where x t is the system microscopic configuration at time t, L(x) is the macroscopic length corresponding to microscopic configuration x, W t is the thermodynamic work done on the system by the external potential, up to the time t, defined by W t = t 0 dt ′ ∂H/∂t ′ , and the average is over all the trajectories of fixed duration t. In order to recover the partition function Z(L) from eq. (5) we use the procedure introduced and discussed in ref. [8]. The investigation of the free energy landscape of the 1TIT molecule, by using Eq. (5), turned out to be a very difficult task. Indeed, the typical value of the work associated to the unfolding of this molecule is of the order of some hundreds of k B T , for the value of ǫ here used. Since one has to evaluate exp(−βW t ), in order to exploit Eq. (5), the JE cannot give a reliable estimate of the energy landscape of the 1TIT molecule. For a discussion of the range of applicability of the JE to microscopic systems see, e.g., [25]. Therefore we consider the PIN1, a smaller protein whose folding characteristics have already been studied with the WSME model [18] (pdb code 1I6C, 39 aminoacids, ǫ/k B = 44 K [21]). Its free energy landscape as a function of the molecule elongation L, as given by Eq. (5), is plotted in Fig. 4 for different velocities of the pulling protocol. It can be seen that the curves F(L) collapse onto the same curve, as the pulling velocity is decreased. This is a clear signature that the energy landscape is correctly reconstructed, and its best estimate is the collapse curve, as discussed in Ref. [10]. On the other hand, the model introduced here is simple enough to allow the exact computation of the function Z(L), and hence we can obtain the exact value of the function F(L): the agreement with the landscape evaluated by the pulling manipulations is found to be very good, see Fig. 4.
In conclusion, we have introduced and studied a model of proteins under external loading. The unfolding length of the titin model is found to be in good agreement with the experimental one. We believe that this result represents a remarkable validation of the model that we have introduced: it suggests that our model, although minimal, captures the basic mechanisms underlying the unfolding process of proteins.
For a small protein, the extended form of the Jarzynski equality gives an estimate of the free energy landscape which is in good agreement with the expected one. We believe that our model can be successfully used to study the interplay between the protein structures and the kinetics of unfolding and refolding under external loading. As an example, by using computer simulations and applying an external force, one can easily determine which are the main contacts stabilizing the molecular structures. Furthermore, our model is also suitable to study the thermal unfolding of the proteins, in order to compare the thermal and the mechanical unfolding paths.
Appendix to "An Ising-like model for protein mechanical unfolding" In this appendix we briefly review the definition of the parameters ǫ, ǫ i j , ∆ i j in the WSME model, and discuss the definition of the new parameters l i j which have been introduced in the main text.
The parameters ∆ i j and ǫ i j , appearing in eqs. (1)-(3) of the main text are chosen following Ref. [1], starting from the protein native structure, as given in the Protein Data Bank (pdb in the following, http://www.pdb.org/). An atomic contact is present (∆ i j = 1) if, in the native state of the protein, at least two atoms from residues i and j + 1 (with j + 1 > i + 2) are closer than 4 Å. In this case ǫ i j is taken to be equal to kǫ, where k is an integer such that 5(k − 1) < n at ≤ 5k, and n at is the number of atomic contacts. As an example, in table I, the values of the contact parameter k, defined as k = ǫ i j /ǫ are listed for the 1TIT molecule.
The lengths l i j , in a generic N + 1 aminoacid protein, are defined as follows. Let us represent the aminoacid i with its N i − C α,i − C i sequence. Taking the native state as the reference configuration, l i j is chosen as the distance between the midpoint of the C i and N i+1 atoms and the midpoint of the C j and N j+1 atoms, see figure 5. If j = i + 1, l i,i+i corresponds to the length of a single aminoacid. If i = 0 the first point is substituted with the position of the N 1 atom, while when j = N + 1 we take the position of the C N+1 atom.
In order to fix the energy scale ǫ, we define the dimensionless unfolding temperature of our model T m , at zero force, as that temperature where the molecule order parameter m = 1/N k m k is equal to 1/2. For a given molecule, using the experimental value of the melting temperature T m , defined as the temperature where half of the sample is unfolded, we define the energy scale as ǫ ≡ T m k B / T m . This approach was used in Refs. [2,3].
Here and in the main text we indicate the proteins either with their common name or with their pdb code.
We find that the midpoint temperature of the 1TIT model molecule, at zero force, takes the value T ≃ 8.03. The experimentally measured unfolding temperature for such a protein is T m ≃ 346 K [4]. Consequently, we choose the energy scale to be ǫ/k B = 43.1 K.
The midpoint temperature of the PIN1 molecule (pdb code 1I6C), at zero force, is found to be T ≃ 7.55. The experimentally measured unfolding temperature for such a protein is T m ≃ 332 K [5]. Consequently, we choose the energy scale to be ǫ/k B = 44 K.