format: text

text:
# Recurrent oligomers in proteins – an optimal scheme reconciling accurate and concise backbone representations in automated folding and design studies

## I Introduction

One of the most fundamental and still unsolved problems in biology is the elucidation of the folding process, that is how a protein sequence undergoes the structural rearrangements that eventually lead to the biologically active conformation (believed to be the free energy minimum) . Since the early studies of Levinthal it was clear that the dynamics of folding to the native state could not be governed by mere random processes ; indeed modern folding theories explain fast folding processes by invoking nucleation-condensation mechanisms or funnel-like energy landscapes  that dramatically reduce the space of visited conformations . Topologic, steric and chemical features are so effective in reducing the space of viable conformations that even at a local level, only few degrees of freedom per amino acid are observed. This fact, originally observed by Ramachandran  has been lately used in a variety of numerical schemes. In these approaches proteins are modeled as chains of one or two interacting centers (representing individual amino acids) with a limited set of local degrees of freedom, such as torsion angles or Cartesian positions, chosen to provide optimal compromises between accurate representation and number of degrees of freedom . These models appear excellent from many points of view with the exception that they fail to capture correlations between torsion angles along the peptide chain.

In this paper we address this problem and propose an optimal way to extend the original idea of Ramachandran of limiting the degrees of freedom of individual residues to strings of consecutive amino acids, showing that they are far from independent. Indeed, their correlations are so strong that, as originally pointed out in a paper by Alwyn Jones and Thirup , it is possible to construct a small data bank of protein fragments that can be used as elementary building blocks to reconstruct virtually all native protein structures. We start by following the seminal idea of Unger et al.  that oligomers of a given length found in a coarse grained representation (such as $`C_\alpha `$ coordinates) of native structures do not vary continuously but they gather in few clusters. Each of these can be represented by a single element (that we term “oligon”) that optimally catches the geometrical and topological properties of the entire basin.

Our approach differs from previous work on the classification of structural fragments  in that the procedure we follow to select the oligons has been explicitly optimised for use in fully automated contexts, especially folding and design attempts . Indeed, such attempts are commonly framed within numerical problems of minimizing suitable functionals (such as energy scoring functions) in structure space. The addition of local constraints reduces drastically the space of viable structures and is undoubtedly a desired feature allowing to keep to a minimum the side-effects of using imperfect parametrizations of the free energy or imperfectly known interaction potentials .

The selection strategy we propose is free of subjective inputs or biases and exploits the full knowledge-based information intrinsic in our data-bank of non-redundant protein structures. An appealing feature of the suggested method is that representative fragments are singled out in order of importance, that is according to the frequency in which they appear in natural proteins. We carry out a series of thorough checks and validations of the clustering strategy and show that the optimal sets of oligons do not suffer from finite-size effects of the data bank. It is shown that the optimal representatives have length equal to 5 or 6 and that with only a few tens of them it is possible to fit virtually any protein within about 1.0 Å co-ordinate root mean square deviations (cRMS) per amino acid. Some of the ramifications of this study are discussed and outlined through preliminary investigations in Section III. The optimal sets of oligons presented and discussed here are made publicly available at http://www.sissa.it/~michelet/prot/repset.

## II Methods and Results

The first step in the creation of a set of optimal representatives is the set up of a sufficiently large data bank of protein structures. Such data bank should cover as best as possible the variety of distinct protein structures observed in nature. At the same time it is important to eliminate correlations and biases in the data bank resulting, for example, from structural homology  For these reasons we compiled our data bank by choosing 75 single-chain proteins from a carefully compiled list of non-redundant structures .

The proteins, listed in Table I were chosen from the SCOP database of non-redundant single-chain proteins covering the most common families: all $`\alpha `$, all $`\beta `$, $`\alpha \beta `$ and $`\alpha +\beta `$ and the most common chain lengths. This ensures that, a priori, the selected structures represent a broad spectrum of structural instances with the least bias or redundancy. As discussed later, the results confirm a posteriori, that the size and quality of the data-bank was sufficient for all practical purposes. Each of the proteins of Table I was partitioned in all the possible fragments of $`l`$ consecutive residues. We considered values of $`l`$ ranging from 3 to 10, for which there are 10936 to 10411 distinct fragments. As in previous studies involved with structural classifications, we retained only the $`C_\alpha `$ coordinates of each fragment .

This approach, is practical and consistent with the idea of having an optimal, but schematic representation of structures. Moreover it is “reversible” to a great extent since the whole peptide atomic geometry can be recovered from the mere knowledge of $`C_\alpha `$ co-ordinates . In turn, if needed, optimal side-chain rotamer positions could be satisfactorily obtained by exhaustive or stochastic methods .

### A Theory: clustering algorithm

The goal pursued in this work is to provide a synthetic, but exhaustive, classification of inequivalent local structural motifs to be used in contexts where a broad exploration of the space of viable protein structures is concerned. Hence, the approach pursued here differs from studies aimed at selecting a restricted number of motif classes to be used in homology modelling or automated recognition/classification of secondary motifs . This distinct goal is accordingly pursued with a novel strategy for the identification of classes that is reminiscent of the clustering technique used by Lacey and Cole in an unrelated context . In the following we shall try to propose a strategy able to perform an optimal subdivision into classes of similarity and, for each of these, provide the best representative. Two of the points of force of such method are the absence of any subjectivity or human supervision through the extensive use of optimal knowledge-based classification criteria and also the fact that similarity classes are automatically extracted and ranked according to their frequency of appearance in natural proteins. This wealth of knowledge-based information provided by the procedure allows to choose the representative set that best matches one’s needs. The clustering procedure we used to partition the fragments in suitable similarity classes is conveniently illustrated by the two-dimensional example of Fig. 1 where 1000 points have been assigned randomly to 4 distinct clusters with the same radius but different size (i.e. number of members). Considerable information about the clusters can be obtained by analyzing the histogram of the distance between all pairs of members in the set. At the simplest level the histogram analysis can reveal two distinct scenarios: a) no clusters are present or b) there are clusters with comparable size and degree of internal similarity. In the first case the histogram distribution is expected to crowd around an average value in a bell-shaped fashion. In the second case, two distinct peaks should occur: one corresponding to the typical distance within classes, the other centered around the (larger) average distance of pairs of members from distinct classes. In the case of very few \[many\] classes, the first \[second\] peak dominates.

The inset of Fig. 1 shows the pair-distance histogram for the set of points in our example. It is evident that it consists of two peaks: the first one extends till about the radius of clusters while the location of the second peak coincides with the typical cluster-cluster distance.

Our goal is to exploit the information obtained from the pair-distance histogram to identify first how many different clusters there are and secondly the optimal representative of each cluster. To do so we follow the intuitive expectation that the best representative of each cluster is the one closest to the cluster center. A deterministic way to identify the center of homogeneous clusters, is to find the member with the largest number of other points within a suitably chosen similarity cutoff (we shall term this number “proximity score”). Indeed, points further from the center will have fewer neighbors. This “election” mechanism is reliable for large and homogeneous clusters.

Hence, we start by choosing the first representative of the set as the one with the highest proximity score. This identifies simultaneously both the largest cluster and its representative. Next, we remove the representative and its cluster from the set and recalculate the proximity score of the remaining points and again we select the member with the highest score. As before we removed it and its cluster and proceed in this iterative fashion until the set of surviving points is exhausted.

When such scheme is applied to the set of Fig. 1 – using a similarity cutoff equal to $`R`$=1 – the optimal representatives of the four clusters (marked with squares) are immediately found and ranked according to their cluster size.

### B Results

We applied the same scheme to analyze our data bank of thousands of protein fragments. This time, the points of the previous example are replaced by the fragments themselves, while the notion of eucledian distance between two points is substituted by the cRMS distance of two fragments , $`X`$ and $`Y`$ of equal length, $`N`$,

$$\sigma (X,Y)=\sqrt{\frac{\underset{k=1}{\overset{N}{}}\left|\stackrel{}{r}_k^{C_\alpha }(X)\stackrel{}{r}_k^{C_\alpha }(Y)\right|^2}{N}}.$$
(1)
This notion of distance is meaningful provided that $`X`$ and $`Y`$ have been previously optimally superimposed with the standard Kabsch procedure  . The calculation of the cRMS of each distinct pair of fragments is the most computationally demanding step since it requires an application of the Kabsch algorithm  for each distinct pair of fragments (e.g. this translates in well over ten million pairs of fragments for lengths of the order of 5).

The histogram of all cRMS of pairs of fragments of lengths in the range $`3l10`$ is given in Fig. 2. It can be seen that all distributions show two distinct peaks, with the exception of $`l=3`$, which appears to be exceptionally short, and hence will be omitted from further analysis.

For the smaller lengths, the first peak collects a substantial amount of “hits”, proving that it is meaningful to assume the presence of classes of similarity. It also appears that the height of the first peak constantly decreases with increasing $`l`$. This confirms the intuition that, by considering very large values of $`l`$ every fragment will be a class for itself. Indeed, for lengths greater than 6, the first peak is hardly discernible from the background. Hence, the mere visual inspection of histogram distributions shows that it would not be justifiable to force the introduction of classes of similarity for lengths above 6. Nevertheless, we shall often present results also for length 7 for the purpose of showing how several unrelated criteria indicate such length as a border-case of viable oligons. An important observation for our subsequent analysis is that the extension of the first peak (the intra-cluster one) depends only weakly on $`l`$ and is about 0.65 Å. This provides an unbiased measure for the similarity cutoff and hence we adopted it. The location of the second “background” peak in the histogram of Fig. 2 gives an estimate of the similarity between unrelated fragments and, hence, crresponds to the cRMS deviation of a random pair of segments. This random pair distance increases with the chain length, but it always well above the value of 2 Å, thus justifying a posteriori the use of similarity cutoffs of the order of 1 Å considered in previous studies .

The advantage of the clustering scheme introduced and used here is that, with modest computational effort (the cRMS distances need to be computed once for all) one has simultaneously both the subdivision in clusters and their optimal representatives. An extra payoff of this approach over other clustering schemes is that the representatives are singled out in order of importance. It is important to stress that there is no stochastic element in the analysis since the assignement of elements to clusters follows a “greedy” deterministic approach. One particular instance where the suggested strategy may fail, is when the “fringes” of distinct clusters overlap, that is when an element falls in the similarity basin of more than one representative. In this situation, more sophisticated clustering techniques (such as those based on k-means analysis ) ought to be adopted in place of the present one, in fact, the iterative removal of assigned members would affect both the choice of the representative and also its score. Although we cannot rule out the presence of fringe overlaps in our data bank, we can exclude it has any substantial significance. Indeed, we have checked that the typical cRMS of the extracted representatives matches the random pair distance which, being much greater than 0.65 Å, makes overlaps highly improbable.

Our analysis identified only 28 representatives for length $`l=4`$, 202 for $`l=5`$, 932 for $`l=6`$ and 2561 for $`l=7`$. As we mentioned before, the existence of a limited repertoire of local folds is a consequence of the existence of a discrete number of degrees of freedom per amino acids, as pointed out by the seminal studies of Levitt on $`n`$-state models . The results obtained here contain significantly more knowledge-based information since, for instance, they also yield the representation score of each representative. It appears that the representation weight (i.e. proximity score) of the fragments decreases very rapidly with the rank (see Fig. 3 and Table II). This is an extremely important feature since it indicates that one might discard the representatives with negligible score and hence work with a subset of the whole data bank. This issue is examined in the next section.

One may expect that the best representatives should belong to the most common structural motifs such as helices, strands or turns, whereas the less frequent ones should correspond to the atypical parts of proteins (structure exceptions). This expectation is confirmed by inspection of the actual shape of the highest ranking fragments; the consensus with the work of  and Unger  shows the reliability with which the main motifs can be identified in different contexts or with different methods. The first four oligons for $`l=5`$ and $`l=6`$ are shown in Figs. 4 and 5 (the structural data for the complete sets is available at the URL given in the introduction). The native environment of the first ten fragments of length 5 and 6 are given in Table II. For $`l=5`$ we have also shown the native environments of the best representatives in Fig. 6. A striking outcome of the clustering analysis is that the first 15 oligons of length 5 and 6 represent over 75 % and 47 %, respectively, of the whole data-bank fragments! To the best of our knowledge, this are the smallest sets of representative fragments able to cover most local structural instances with an uncertainty comparable with the best experimental resolution.

## III Discussion

### A Analysis of the clustering procedure

Before testing the goodness of the representative fragments it is necessary to validate the clustering procedure and ensure that the results are robust and not too dependent on the details of the data bank. We carried out a first check by studying how the outcome of the clustering scheme is affected by the size of our data bank. To be precise, this test goes beyond the mere validation of the oligon extraction scheme, since it also constitutes a check of the applicability of any clustering scheme to protein fragments. To proceed in an unbiased way we randomized the order of fragments in the data bank, so to cancel correlations of consecutive (overlapping) oligons, and extracted the representatives for an increasing number of fragments taken from the top of the randomized list. A careful analysis of the data has revealed that, for any length, $`l`$, the number of trivial representatives, i.e. those that, having score equal to 1, represent only themselves, grows linearly with the size of the data bank. The preportion of trivial representatives is about 0.5 %, 2% of the whole population for lengths 5, 6. For length greater than 6 the proportion of trivial representatives is considerable (being greater than 10 %). On the other hand, for $`4l6`$ the number of non-trivial representatives shows very little increase with the size of the data bank and can be considered constant for all practical purposes. This provides a solid a posteriori confirmation that the data-bank is of sufficiently large size. Of course, the number of representatives and their growth with data-bank size depends on the particular choice of similarity cutoff (the smaller the value, the larger the number of classes). In this particular study the choice of the cutoff was dictated by the properties of the very same data to be clustered. Nevertheless, the use of physically viable cutoffs lead, invariably, to the identification of the same high-ranking clusters and, correspondingly, almost identical representatives. This could be expected a priori since identifying the most common local folds should be independent, to a large extent, on the details of the clustering procedure. As explained in the next sectiom, we tried to build on this robust result and concentrate only on the top representatives.

### B Reducing the representative sets

Since the trivial oligons mentioned in Sec. III A represent only themselves, one may wonder if they can be dropped from the set and still be able to represent well the majority of native structures. In this subsection we considered this problem and try to quantify the attainable accuracy in representation when a subset of the representatives is used. For a preassigned number, $`m`$, of representatives to be used, the optimal accuracy is obtained when the $`m`$ highest ranking fragments are taken. Hence, our extraction scheme is particularly convenient for this type of study since it yields the representative fragment ranked according to their proximity score

In this framework, we measure the accuracy of representation by the amount of local structural deformation required to bring all fragments of the native structure within the proximity basin of any of the reduced oligons (while in Sect. IIID we shall fit several proteins with the oligons). This is a sort of measure of the “completeness” of the set of oligons: if the used set of oligons represented all possible instances of protein fragments, no deformation would be required. On the contrary, the poorer the set of representatives, the larger is the deformation required to bring the original fragments in the proximity basin of one oligon. To do so, we use a stochastic Monte Carlo dynamics on the backbone (described in Appendix A) to minimise the following quantity:

$$S\underset{i=0}{\overset{L/l1}{}}(\sigma (B_i,\omega _i)R)^2\theta [\sigma (B_i,\omega _i)R]$$
(2)
where $`L`$ is the length of the protein, $`\sigma `$ is the cRMS distance of eqn. 1, $`B_i`$ is the $`i`$th backbone fragment of length $`l`$, $`\omega `$ is its closest oligon, R is the similarity distance (0.65 Å) and $`\theta `$ is the usual step function.

By using the stochastic dynamics, the starting structure is deformed until all fragments are within the preassigned distance R from one of the oligons. When this happens, the score function, $`S`$, is exactly zero and the dynamics is stopped. By measuring how far (in terms of cRMS) the backbone has moved from its original position we can judge whether the achievable quality of representation is acceptable. We carried out this scheme by using only the first few representatives (for each length $`4l7`$ ) and then increased their number progressively (always choosing the highest ranking ones). The cRMS as a function of the number of representatives is shown in Figs. 7. For $`4l6`$ only a fraction of the collected oligons are necessary to fit the 10 test structures within 0.65 Å and with no need to distort them. Even for $`l=6`$ with only 100 representatives any protein backbone can be fitted at the price of minute distortions (less than 0.5 Å ), that are finer than the typical experimental structural resolution.

It is important to point out the the low global cRMS values given in Fig. 7 do not hide exceedingly large local distortions averaged with many other smaller local deviations. Indeed, the deviations appear to be homogeneous along the chain; the worst local distance of fitted fragments from the native positions never exceeds twice the global averaged value (data not shown). To be more precise in assessing the presence and effects of unphysical local deviations we calculated the displacementes of $`C_\beta `$ positions in fitted backbones from the native one. Indeed, $`C_\beta `$’s discrepancies are good indicators of local variations in the dihedral angles between virtual $`C_\alpha `$ bonds. $`C_\beta `$ positions were recovered from $`C_\alpha `$ coordinates $`\{\stackrel{}{r}_i^\alpha \}`$ through the standard geometric constrained construction :

$$\stackrel{}{r}_i^\beta =\stackrel{}{r}_i^\alpha +d_0\left(\widehat{a}\mathrm{cos}(\theta )+\widehat{b}\mathrm{sin}(\theta )\right)$$
(3)
where:

$$\widehat{a}=\frac{\widehat{s}_{i,i1}+\widehat{s}_{i,i+1}}{\left|\widehat{s}_{i,i1}+\widehat{s}_{i,i+1}\right|}\widehat{b}=\frac{\widehat{s}_{i,i1}\widehat{s}_{i,i+1}}{\left|\widehat{s}_{i,i1}\widehat{s}_{i,i+1}\right|}$$
(4)
and:

$$\widehat{s}_{i,j}=\stackrel{}{r}_i^\alpha \stackrel{}{r}_j^\alpha .$$
(5)
In the previous formulae $`d_0=3`$ Å is the distance of the $`C_\beta `$ atoms from the corresponding $`C_\alpha `$ atom and $`\theta `$ is the out-of-plane angle optimally set to $`37.6^0`$. The $`C_\beta `$ positions are very sensitive to the local position of the $`C_\alpha `$ because a wrong (even by a small amount) choice of the angle between the $`\widehat{s}_i`$ can heavily affect its position, e.g. shift it to the wrong side of the chain.

We considered some of the proteins in the test set previously fitted with a subset of the representative oligons. For these we constructed the $`C_\beta `$ positions and calculated the deviations of the latter from those in the native configurations. the data are shown with dotted lines in Fig. 7 and highlight how the discrepancy is very small and follows the trend of the cRMS for $`C_\alpha `$ atoms. This shows that the local distortions are really tiny even when 100 of the over 600 oligons of length 6 are used. As usual, an atypical behaviour is seen for length 7, for which, even using hundreds of fragments, a much larger discrepancy is observable.

### C Optimal length of representative oligons

Each of the sets of representative fragments of length $`4l6`$ are optimal by construction and all of them satisfy the rigid tests carried out so far. The goal we pose here is to decide which length is the best. The answer is certainly not unique, since different criteria for optimality can be used . For example, if one is interested in having the smallest possible set of representatives, then small values of $`l`$ are to be preferred. On the other hand, if one is mainly interested in having the least number of conformational degrees of freedom per residue then $`l`$ should be chosen as large as possible. Both approaches can be legitimate in appropriate contexts. From a general point of view, however, using very short fragments defeats the purpose of this study - that is to capture structural correlations. On the other hand, excessively large values of $`l`$ are more difficult to handle and uninteresting since clusters will typically be sparsely populated (over specialised case). Here we examine the main properties of representative oligons that can be conveniently exploited in different contexts. We begin by discussing how well oligons of different length represent secondary motifs . The latter are indeed the distinctive feature of proteins (as opposed to random heteropolymers  and have several consequences on biophysical properties, such as speeding up the folding process or providing maximum kinetic accessibility to the native state .

Alpha-helices seem to be fairly easy to represent. In fact, for all cases $`l=4,5,6,7`$ a single representative (namely the highest-scoring one) is sufficient to represent virtually all instances of helices. The situation is different for $`\beta `$-strands, due to the different environment in which they can be found (parallel or anti-parallel, bent $`\beta `$-barrels, Greek-key motifs etc). This variability implies that more than one representative for $`\beta `$ motifs is found (although not with the same proximity score) . Examples are shown in Figs. 4 and 5. This proliferation effect is more dramatic for longer fragments, consistently with the findings of Prestrelsky et al. . Indeed, for $`l=7`$, each of the distinct $`\beta `$ classes appears severely depleted, containing typically less than 100 elements which is a small fraction of the score of the helical one (2070).

For all values of $`l`$, however, the largest number of representatives is covered by segments representing loop regions. These results are particularly relevant for modeling/characterizing regions of high variability, but our main focus is on the possibility to represent synthetically, though accurately, recurrent oligons. Within such minimalistic approaches, the choice of representatives of length 5 seems to be the best one, since it captures non-trivial correlations while using essentially a single representative for $`\alpha `$ and $`\beta `$ instances.

### D Fitting proteins with oligons

Another criterion for selecting the most suitable length is how well can we reproduce a given protein by “gluing” rigidly together only the representative oligons? The purpose of such question is to investigate the benefit of employing oligons in folding contexts. A simple and powerful way to speed up the numerical simulations of folding would be to consider structures made only by “gluing” suitably chosen representative oligons. In such framework the only degrees of freedom that one has to contend with are: 1) which oligon to use and 2) how to connect successive oligons. This is a severe reduction of the traditional continuous/discrete degrees of freedom per amino acid adopted in ordinary Monte Carlo or Molecular Dynamics schemes. The feasibility of such scheme depends first of all on the possibility to reproduce sufficiently well any given native structure by joining rigidly the oligons. We checked this by following a stochastic process to find both the best oligons to be used locally and also their best relative orientations. This was almost a worst-case scenario due to the independence of the test set from that of Table I. The optimal fit was accomplished by progressively distorting the native structure with the local Monte Carlo moves described in Appendix A. The ”energy-like” cost function had the same form of (2), but where $`R`$ is set to an arbitrary small positive quantity, $`10^3`$ in our case. Again, we carried out the stochastic dynamics (proceedings through very tiny local deformations) until the cost function was reduced to zero. This signalling that each protein fragment had been optimally collapsed on an oligon. It can be anticipated that, due to the propagation of misfits, the cRMS with respect to the native protein would be rather larger than the similarity cutoff of 0.65 Å. Moreover, it may be expected that smaller oligons may lead to smaller cRMS since they might provide more flexibility in “tiling” target structures. Surprisingly, this is not the case, as visible in Table IV, where we summarised the global cRMS deviations for rigidly fitting the 10 proteins in the test set. Remarkably, the overall cRMS is always very close to 1 Å such cRMS deviations of the native and fitted protein can be appreciated visually in Fig. 8. We explain the little dependence of cRMS fits on oligon lengths with the observation that, irrespective of the oligon length, each residue in native conformations is typically 0.5 Å away from the corresponding position in the best-matching oligon. This little sensitivity on $`l`$ is, in turn, reflected on the overall cRMS of the rigid fit. The fit discrepancy is not only independent of the length but also fully compatible with state-of-the-art experimental resolution of crystallographic structures. For these reasons one may adopt oligons of the longest possible lengths if the primary interest is capturing the longest possible structural correlations. This would suggest to consider lengths equal to 6. Our fit scheme has considerable advantages over previous ones where representatives obtained with different techniques were employed. For example, in their classic paper, Unger et al.  used a molecular best fit procedure that yielded cRMS of over 7 Å when hexamers were used to fit peptides of over 70 residues. The dramatic improvement of the results in Table IV confirms the validity and reliability of both the clustering method and of the extracted set of oligons. Indeed, the low values of cRMS fit support the expectation that the extracted oligons can be successfully used to speed up folding attempts. Preliminary tests in this direction have been carried out in folding contexts where perfectly smooth folding funnels  lead to known crystallographic structures. Such studies originally undertaken to elucidate global aspects of the folding process have recently been the key to predict and describe the influence of topological protein properties on folding nuclei and/or thermodynamical folding stages . By employing oligons of length 5 we were able to speed up the collection of folding data by several factors .

### E Correlation between oligons and amino acid sequences

We devote the final part of this section to elucidate the possibility of finding correlations between oligons and amino-acids sequences. In general, it is well-known that there is preference for definite sets of amino acids to occupy or avoid specific structural motifs . Here we examine the extent to which such propensities are reflected in the oligons and the clusters they represent. Highlighting connections between sequences and oligons has a twofold purpose: a clear preference of an amino-acid sequence to be mounted in a specific oligon can be useful exploited in folding predictions, whereas design attempts can be greatly aided by discovering that some oligons preferably house very few sequences.

The connection between sequence-structures connections have been heavily investigated, with fair success, for a variety of fragment lengths and amino-acid sequences. It is important to examine the issue also in the present context since the emergence of clear correlations between sequences and oligons could be an additional aid in reducing the computational complexity of folding and/or design.

For sake of simplicity we consider in this section only the case $`l=5`$ and we considered the best $`40`$ oligons of that length. We start by introducing a suitable classification of the 20 types of amino acids. This is essential to proceed, since otherwise the shear number of the possible sequences, $`20^5`$ 3 million, would make it impossible to gather sufficient statistics for all quintuplets. The classification scheme we introduce here is based on some general results for chemical affinities  and some empiric attempts. According to it we subdivide the residues in four distinct classes.

In the first we place Gly, in the second Pro, in the third the hydrophobic (H) aminoacids (Ala, Val, Leu, Ile, Cys, Met, Phe, Tyr, Trp) and finally in the fourth the polar (P) ones (Hys, Ser, Thr, Lys, Arg, Asp, Asn, Gln, Glu). With this subdivision we keep separate the amino acids (Gly and Pro) that can attain atypical conformations/chiralities  (and hence may act as helix breakers etc.). It is also wise to keep in separate families hydrophobic and polar aminoacids, since they can alternate regularly in secondary motifs partly exposed to the solvent . Within this framework we could obtain in principle up to $`4^5=1024`$ distinct pentamer sequences (we always consider our pentamers as “directed” in that the C and N termini are not exchangeable). It turns out that, due to chemical and steric constraint, not all pentamer sequences are observed in nature, and hence in our data-bank.

To perform our analysis we considered all the proteins (75) appearing in Table I . We partitioned then in overlapping fragments of length 5 ending up with 10786 pentamers. The size of this data-bank was sufficient to provide excellent coverage of all possible pentamer sequences. This is evident from the plot of Fig. 9 which shows how the number of distinct pentamer sequences grows with the data-bank size.

The asymptotic number of distinct sequences we obtained from the near six thousands instances was 614, about half of all possible ones.

To match the 614 sequences to the 40 oligons of length 5 we re-applied the clustering procedure: to each of the oligons we assign not only its native sequence but also those of each member in the cluster it represents. All this information can be conveniently stored in a score matrix $`z(i,j)`$ whose entries correspond to the number of times that the $`j`$th sequence has been assigned to the $`i`$th oligon (hence $`z`$ is a 40x614 matrix);

A two-dimensional representation of the score matrix is plotted in Fig. 10 where the dark boxes correspond to entries above 25, the grey ones to entries between 3 and 25 and the blank ones to entries below 3 . The figure shows that $`z`$ is a sparse matrix, since only few entries have a significative entry (bigger than 3). This supporting the conjecture of strong correlations between oligons and sequences.

The last observation can be turned into a more quantitative statement by examining the behaviour of definite oligons and/or pentamer sequences. The natural candidates to focus on are the 135 sequences that appear more that 20 times, and hence allow a statistically sound analysis. For each one of these sequences we examined the relative frequency with which they occupy a given oligon. Typical results are given as histograms in Fig. 11.

It appears that sequences do not occupy many oligons; in fact, less than 18 oligons are occupied, on average, by the 135 sequences ( and over 70 % of the entries is covered by six oligons). It is worth underlining how this is not an average effect reflecting the relative magnitude of the proximity scores of the oligons. To show this one can establish a reference threshold corresponding to the number of expected hits if sequences are distributed uniformly over all fragments. Thus, for a given oligon, the threshold is simply the ratio between its proximity score and the total number of fragments used to calculate this score. It was found that in 103 cases out of 135 ( 77 % ) the sequences select their preferred oligon with a percentage significantly higher (in excess of $`20\%`$ than the trivial threshold). Although it is clear that any given sequence is compatible only with few oligons, the converse is not true. This interesting asymmetry between sequence and structure has deep roots, as first shown by Anfinsen , who pointed out that a protein sequence uniquely identifies its structure, while several different sequences can admit (almost) the same structure as their native states. This aspect is strikingly evident when plots analogous to those of Fig. 11 (by interchanging the role of sequences and oligons) are made. In Fig. 12 the occurrence frequency histograms for the first (ranked according to the proximity score) four oligons are plotted. In these histograms for each sequence (listed in ordinate according to a convenient scheme) the percentage of occurrence for the given oligon is represented.

It is clear that, unlike the case for pentamer sequences, there is not a preference for a given oligon to be occupied by few sequences, so that the benefits of these correlations studies for design schemes is not as dramatic as could be for folding simulations.

As a final test we verify whether it is possible to define selection rules for locating amino acids in well-defined oligon positions, e.g. to pinpoint particular points where it is unlikely that some class of amino acid could appear. The existence of such forbidden points could be, again, a useful source of information for folding and design. Due to the non-homogeneous population of the amino acids classes we adopted, we expect to extract information only for the first two classes, namely Gly and Pro. For any oligon we considered all the related sequences and we monitored, site by site, the occurrence frequency of each class. If for a given site and class this frequency is below the threshold of 0.5% we consider the event unprobable (and hence significant in the present context). In Table V we list 19 of these events which take place in the highest-ranking oligons. The most significant instances all refer to Pro “class”. The final test we have summarised shows that a combination of the structural reduction in oligons and associated correlations with local sequence propensities can be turned into a powerful tool in aiding folding and design. This hope is corroborated by the recent successes of structural prediction schemes based on local sequence propensities .

## IV Conclusions

The starting point of this work was the conclusion of recent previous studies that there exist recurrent local motifs in natural proteins . We introduce novel and fully automated criteria for an optimal partitioning of a complete data-bank of fragments taken from non-redundant proteins into classes of similarity.

We exploit the intrinsic information in the data-bank to identify the classes with the least bias or human supervision. Our goal was to show that such scheme succeeds in conciling two competing aspects of protein modeling: accuracy and synthetic modeling .

In fact, on one hand this method is shown to provide the most economic subdivision in classes (the number of which is not set a priori). On the other hand, the optimally extracted representatives from each class are shown to be sufficient to represent and fit virtually all protein structures with an uncertainty of 1 Å(rigid fitting) or 0.5 Å, when only local similarity within the proximity basin is required. We also considered several possible lengths for oligons and examined their suitability in different modelling contexts. It turns out that $`l=5`$ is the most suitable when the smallest representative set is needed, while $`l=6`$ is best when it is necessary to capture the longest possible correlations. Lengths smaller than 5 or longer than 7 appear to be far from optimality.

## V Acknowledgments

We acknowledge support from the Theoretical and Biophysical sections of the INFM. We are indebt to the Italian Research Council for the financial support of the advanced research project “Statistical mechanics of proteins and random heteropolymers”.

## A Monte Carlo dynamics

In this Appendix we present a summary of the stochastic approach that we used for the dynamics of the protein backbones. As mentioned in the text we used the Monte Carlo dynamics for progressive distortion of native protein backbones in order to fit them locally by using a restricted set of representative fragments (see Section III) to provide the best protein fit by using exactly the representatives. In the spirit of standard dynamical approaches for three-dimensional structures  each time we propose a Monte Carlo move we distort the structure by performing either local or global rearrangements. Local moves are single-bead or crankshaft, as explained below, while pivot rotations were employed for global ones.

In the following we will use the ordinary Cartesian triplet $`(x,y,z)`$ to indicate the co-ordinates of $`C_\alpha `$ atoms. Subscripts will denote the amino acid position along the sequence. The three types of moves are as follows:

1. Single $`C_\alpha `$ move.  A random site $`i`$ of the protein chain is chosen and its old coordinates are replaced by new ones $`(x_i^{^{}},y_i^{^{}},z_i^{^{}})`$ defined as:

$$x_i^{^{}}=x_i+\eta _1\mathrm{\Delta }ly_i^{^{}}=y_i+\eta _2\mathrm{\Delta }lz_i^{^{}}=z_i+\eta _3\mathrm{\Delta }l$$
(A1)
where $`(\eta _1,\eta _2,\eta _3)`$ are three independent random numbers in the interval $`(1,1)`$ and $`\mathrm{\Delta }l`$ is a distance that we fixed (see discussion below) equal to 1 Å (top panel of Fig. 13).
2. Crankshaft move. Two protein sites $`i`$ and $`j`$ with sequence separation at most $`10`$ are chosen. Then the all the sites between $`i`$ and $`j`$ are rotated around the axis going through $`i`$ and $`j`$ by a random angle in the range $`\frac{\pi }{10}\theta \frac{\pi }{10}`$; (middle panel of Fig. 13).
3. Pivot move. A random site $`i`$ and a random axis passing through it are chosen. All the sites from $`i+1`$ to the end are then rotated around the axis by an angle in the range $`(\frac{\pi }{10},\frac{\pi }{10})`$; (bottom panel of Fig. 13).

The new configuration generated by applying one of these moves (chosen with equal weight) is first examined to make sure that it does not violate basic geometrical constraints obeyed by natural proteins, namely:

1. the distance between two consecutive $`C_\alpha `$ atoms (measured in Å ) must remain in the range $`(3.7,3.9)`$ and
2. the distance between two non consecutive $`C_\alpha `$ atoms must be greater than $`4\AA `$.

If these conditions are not fulfilled, then a new move is attempted. When the new configuration has passed the geometrical test then is accepted/rejected through the classic Metropolis rule.


text:
# Complex spatial organization in a simple model of resource allocation

## Abstract

A dynamical model for the distribution of resources between competing agents is studied. While global competition leads to the accumulation of all the resources by a single agent, local competition allows for a wider resource distribution. Multiplicative processes give rise to almost-ordered spatial structures, thourgh the enhancement of random fluctuations.

Though many real systems in the scopes of biology and the social sciences are well described by agent-based models with pair interactions , qualitatively similar to the physical description of interacting particles, a large class of biological and social processes are driven by interactions mediated by an external actor, which generally bears no resemblance with the individual agents. This kind of processes—which, in connection with physics, may be assimilated to the evolution of globally coupled dynamical systems —is typically found in systems where agents compete for resources. A particularly relevant instance of this situation is given by competing biological species. In this case, in fact, the interaction is rarely given by struggling events between individuals of different species, but rather by accessing simultaneously to the limited resources provided by the environment. The performance of each agent, i.e. of each species, is here measured by its ability to get resources, which depends both on the capabilities of each individual and on global features such as the total population of the species. The same scenario is found in some economical systems, for instance, in companies competing for financial resources, usually administrated by banks, or even in scientific research projects competing for funds from a government agency. A key ingredient in the dynamics of these systems is that the ability of each agent to get resources at a given time depends, often strongly, on the resources assigned to the agent at previous stages. This can give rise to a multiplicative process that, in the absence of buffer mechanisms, leads to resource accumulation by a single agent. In the following, we consider a simple dynamical model that incorporates these elements and, in particular, study the effects of competition at local level. We find that local competition softens the process of resource accumulation and gives rise, through the enhancement of spatial fluctuations, to nontrivial structures.

Our system consists of an ensemble of $`N`$ agents, each of them characterized by its productivity $`\alpha _i`$. For future convenience, we consider the generic situation where all productivities are different. At each time step $`t`$, each agent is assigned an amount of resources $`r_i(t)`$, which is used to produce an amount of products $`p_i(t)=\alpha _ir_i(t)`$. At the next time step, resources are distributed among agents in amounts proportional to their production, namely,

$$r_i(t+1)=\frac{p_i(t)}{_jp_j(t)}R(t+1)=\frac{\alpha _ir_i(t)}{_j\alpha _jr_j(t)}R(t+1),$$
(1)
where $`R(t+1)`$ are the total resources available at time $`t+1`$. Model (1) can be fully solved for arbitrary productivities and initial conditions $`r_i(0)`$. In the first place, we note that rescaling $`r_i(t)/R(t)r_i(t)`$, Eq. (1) becomes independent of the total resources. We thus fix $`R(t)=1`$ for all $`t`$, so that $`_jr_j(t)=1`$. In this situation, the solution to Eq. (1) reads

$$r_i(t)=\frac{\alpha _i^tr_i(0)}{_j\alpha _j^tr_j(0)}.$$
(2)
For asymptotically large times, $`r_i(t)0`$ for all $`i`$, except for the agent with the maximal productivity, $`\alpha _{\mathrm{max}}=\mathrm{max}_i\{\alpha _i\}`$, which receives all the available resources. Due to the multiplicative effect of resource allocation according to production, all resources are in the long run accumulated by the agent with the maximal productivity, giving rise to a sort of winner-takes-all state . By analogy with population dynamics we say that the remaining agents become extinct. In fact, in connection with biological populations, this result is a realization of a well-known principle of ecology, namely, the principle of competitive exclusion : in a system of biological species competing for the same resources, only one survives and the others undergo extinction. At moderately large $`t`$ the resources assigned to each agent are well approximated by $`r_i(t)=(\alpha _i/\alpha _{\mathrm{max}})^tr_i(0)/r_{\mathrm{max}}(0)`$. Assuming that the productivities $`\alpha _i`$ are drawn at random from a distribution $`P_\alpha (\alpha )`$ and that, for simplicity, resources are evenly assigned at the beginning, $`r_i(0)=N^1`$ for all $`i`$, the probability distribution for the individual resources at a given time is

$$P_r(r_i)=\frac{\alpha _{\mathrm{max}}}{r_it}r_i^{1/t}P_\alpha \left(\alpha _{\mathrm{max}}r_i^{1/t}\right).$$
(3)
The dependence of this function on $`r_i`$ through the power $`r_i^{1/t}`$ is very weak for large $`t`$, so that on a wide interval of the variable we find $`P_r(r_i)r_i^1`$. For long times, thus, resources are distributed among agents following a power law with exponent $`1`$. Compare this result with Pareto’s law of wealth distribution .

Model (1) admits several variations, that may be of interest in connection with the description of some real systems. For instance, the extinction of all but one agent can be avoided by modifying slightly the allocation method. If, at each time step, a small fraction $`\rho `$ of the total resources is evenly distributed among the agents whereas the remaining is assigned according to production as above, the individual resources are always larger than $`\rho /N`$. In this situation, $`P_r(r_i)`$ becomes stationary for long times and, if $`\rho `$ is small enough, the power-law dependence quoted above is found again for $`r_i>\rho /N`$.

A second variation consists in assuming that the productivity can in turn depend on the individual resources, for instance, as $`\alpha _i(r_i)=\alpha _i^0A(r_i)`$. If $`A(r)`$ increases with $`r`$, the winner-takes-all effect is enhanced. The opposite case, where $`A(r)`$ is a decreasing function, is more interesting. Lower productivities for higher resources may be the consequence of size effects, crowding, or loss of efficiency due to lower competition. In this case, the system evolves to a stationary situation where the agents whose productivity weights $`\alpha _i^0`$ are above a certain threshold $`\alpha _{\mathrm{min}}^0`$ receive nonzero resources, whereas the agents with $`\alpha _i^0<\alpha _{\mathrm{min}}^0`$ become extinct. In the stationary situation, all the surviving agents have the same productivity $`\alpha `$. The third instance, where $`A(r)`$ is a nonmonotonous function, can give origin to multiple stationary nontrivial solutions.

In the following we focus the attention on a variation of model (1) that introduces a spatial distribution of agents. Though the total resources are still allocated among all the ensemble, the agents compete at a local level only. The set of agents that compete with a given agent $`i`$ defines its neighborhood $`𝒩_i`$. Resources are assigned according to

$$r_i(t+1)=\frac{1}{Z(t)}\frac{\alpha _ir_i(t)}{_{j𝒩_i}\alpha _jr_j(t)},$$
(4)
where

$$Z(t)=\underset{i}{}\frac{\alpha _ir_i(t)}{_{j𝒩_i}\alpha _jr_j(t)}$$
(5)
is a normalization factor that insures that $`_jr_j(t)=1`$ for all $`t`$. We have performed numerical simulations of this system in ensembles of $`N=10^3`$ to $`10^4`$ agents distributed on various geometries. The productivities $`\alpha _i`$ were drawn at random from a uniform distribution in ($`0,1`$). Note that, in Eqs. (1) and (4), a homogeneous rescaling of all productivities leaves the models invariant. Several distributions of initial conditions were tested, but there are no essential differences with the case where, at the first step, resources are evenly allocated, $`r_i(0)=N^1`$. Therefore, we concentrate our simulations on this simple case.

In the first place, we consider a one-dimensional array of agents with periodic boundary conditions. The neighborhood of the $`i`$th agent consists of its two nearest neighbors. We find that, for long times, approximately half of the agents become extinct and resources are evenly allocated among the surviving agents. Practically everywhere along the array agents are ordered in an alternating sequence of extinct agents and survivors. Occasionally, one finds two neighbor sites where both agents are extinct, but two survivors are never contiguous. In other words, in the neighborhood of a surviving agent no other survivor can be found \[cf. model (1)\]. The surprising feature of this asymptotic distribution is that the strong correlation between survival and productivity found for Eq. (1) is apparently lost in Eq. (4). In fact, productivities are distributed completely at random along the array, whereas the resulting distribution of survivors is highly ordered. At first glance, no connection can be established, for instance, between survivor sites and local maxima of the productivity or any spatial pattern caused by fluctuations in the distribution of $`\alpha _i`$.

The same feature is found in a two-dimensional array with periodic boundary conditions, where the neighborhood of each site consists of the four nearest neighbors. Figure 1 displays the asymptotic state on a $`100\times 100`$-site lattice ($`N=10^4`$), where the full squares represent survivor sites. Regular domains where survivors and extinct agents alternate in both dimensions, separated by worm-like boundaries formed by extinct agents, are apparent. We have also verified that different definitions of the neighborhood of a site, both in one and in two-dimensional lattices, produce similar almost-periodic asymptotic structures. Again, resources are evenly distributed between the surviving agents and, in neighborhood of a survivor, no other survivor can be found. For instance, in a one-dimensional array where the neighborhood consists of four sites (nearest and next-to-nearest neighbors), the resulting structure is a periodic sequence of two extinct agents and one survivor. Occasional defects, with larger zones of extinction, are also found.

In order to trace the origin of the almost-periodic structures emerging at asymptotically long times we have inspected the successive states of the system from the earliest evolution stages, for the specific case of a one-dimensional array of agents with two neighbors per site and periodic boundary conditions. As above, productivities are drawn at random from a uniform distribution and $`r_i(0)=N^1`$ for all $`i`$. We find that even at the first step of the evolution the distribution of resources already exhibits an almost-periodic sequence of relatively high and low values, in spite of the fact that productivities are spatially uncorrelated. Namely, for a substantial fraction of agents, we have $`[r_{i+1}(1)r_i(1)][r_i(1)r_{i1}(1)]<0`$. This observation reveals that the simple dynamics of model (4) is able to introduce strong spatial correlations to an initially uncorrelated state, even at the level of a single evolution step.

To quantify this effect, we take a set of five uncorrelated random numbers $`\alpha _1,\mathrm{},\alpha _5`$ drawn from the productivity distribution, and consider the combinations

$$\sigma _i=\frac{\alpha _i}{\alpha _{i1}+\alpha _i+\alpha _{i+1}}$$
(6)
for $`i=2,3,4`$. The quantities $`\sigma _i`$ are then proportional to the resources received at $`t=1`$ by three consecutive agents in the array \[cf. Eq. (4)\]. Numerical realizations of these quantities, over $`10^6`$ independent choices of $`\alpha _1,\mathrm{},\alpha _5`$, show that $`(\sigma _4\sigma _3)(\sigma _3\sigma _2)<0`$ with probability $`p_1=0.762`$. Consequently, slightly more than 76 % of the agents are expected to belong to a periodic sequence of alternating high and low resources at the first time step. The remaining 24 % stands for defects in the periodic structure. The same kind of analysis can be performed for successive steps in the evolution. It is found that correlations are further enhanced by the dynamics. At the second evolution step, for instance, the probability for an agent to belong to the periodic structure grows to $`p_2=0.782`$. For later stages, we find $`p_5=0.816`$, $`p_{10}=0.839`$, $`p_{20}=0.854`$, and $`p_{50}=0.864`$. These probabilities saturate at $`p_{\mathrm{}}0.865`$.

While the emergence of almost-periodic structures seems to indicate that no correlations subsist in model (1) between survival and productivity, the inspection of the evolution equations as well as of Eq. (6) suggests that some remaining correlation should however exist. In fact, low productivities still imply few resources, so that the minima in the periodic structure should preferably coincide with sites with small $`\alpha _i`$. Moreover, defects in the periodic structures—where extinct agents are found in consecutive sites—should also correspond to low-productivity zones. Such remaining correlation can be characterized by the probability of survival as a function of the productivity, $`p_s(\alpha )`$. This probability is defined as the fraction of agents with productivity $`\alpha `$ which survive at asymptotically long times. We have numerically measured $`p_s(\alpha )`$, in series of $`100`$ independent realizations of the productivity distribution for ensembles with $`N=1000`$, for various geometries and neighborhoods. Figure 2 shows the survival probability as a function of $`\alpha `$ for one-dimensional arrays with different numbers of neighbors. Here, the correlation between survival probability and productivity is apparent. Agents with larger productivities are more likely to survive that those with small $`\alpha `$. Note, however, that for $`\alpha =1`$ the survival probability is less than unity, so that even with the maximal productivity there are chances of undergoing extinction. Conversely, for $`\alpha 0`$ (but $`\alpha 0`$), the survival probability is finite. As the number of neighbors grows the probability is more concentrated towards larger productivities, as expected. In the limit where the neighborhood extends to the whole array the original model (1) is recovered, and the survival probability must vanish except for $`\alpha =1`$.

The survival probability can be calculated independently of the numerical realization of the model, using expressions as in Eq. (6). For a fixed value of $`\alpha _3`$ and different random choices of $`\alpha _1`$, $`\alpha _2`$, $`\alpha _4`$ and $`\alpha _5`$, an estimate at the first evolution step of the survival probability of the agent with productivity $`\alpha _3`$ is given by the fraction of instances for which $`\sigma _3>\sigma _2,\sigma _4`$. Generalizing this procedure for successive time steps, the estimate can be improved by considering later stages in the evolution. The line in Fig. 2 corresponds to this estimate at $`t=50`$. It shows an excellent agreement with the numerical realizations.

To compare now several geometries, we have measured the survival probability in systems where each agent has the same number of neighbors—specifically, four—on different kinds of arrays, namely, a one-dimensional array (nearest and next-to-nearest neighbors), a two-dimensional lattice (nearest neighbors) and a random graph with four connections per site. The results, shown in fig. 3 corresponds to averages over $`100`$ realizations in systems with $`N=1000`$ for one dimension and random graphs. In this latter case, the graph topology is chosen anew at each realization. For two dimensions, the results correspond to $`10`$ realizations on a $`100\times 100`$-array. Though the general trend of $`p_s(\alpha )`$ is qualitatively the same for the three cases, some systematic differences are apparent. For instance, the survival probability in one dimension is appreciably lower than in the other geometries for low productivities. In this range, remarkably, the data for random graphs lies between those for one and two dimensions. For $`\alpha >0.8`$, instead, the values of $`p_s(\alpha )`$ are hardly distinguishable within our numerical precision.

We may summarize our main results as follows. First, the strong, deterministic correlation between survival and maximal productivity that characterizes the model with global competition (1) is replaced, in the model with local competition (4), by a much weaker, probabilistic-like correlation. Agents with high productivity are more likely to survive but, even with the maximal productivity, extinction cannot be completely discarded. On the other hand, agents with very low productivity have a chance of survival. This conclusion should be relevant to the possible applications of the present models, both to economy and to biology. Second, the loss of the above mentioned deterministic correlation is accompanied by the formation of an almost-regular spatial structure, with a periodic alternation of survivors and extinct agents. This structure is explained by the emergence of strong spatial correlations out of the fully uncorrelated distribution of productivities, due to the very action of the evolution rules. Nontrivial correlations in the conditional probabilities for certain combinations of uncorrelated variables—a rather simple but counterintuitive phenomenon—have recently been pointed out for series of random numbers . The present model illustrates the occurrence of the same kind of phenomenon in a spatially extended dynamical system.

This work has been partially carried out at the Abdus Salam International Centre for Theoretical Physics (Trieste, Italy). The author thanks the Centre for hospitality, and G. Abramson for useful comments.


text:
# Exactly Solvable Single Lane Highway Traffic Model With Tollbooths

## Abstract

Tolls are collected on many highways as a means of traffic control and revenue generation. However, the presence of tollbooths on highway surely slows down traffic flow. Here, we investigate how the presence of tollbooths affect the average car speed using a simple-minded single lane deterministic discrete traffic model. More importantly, the model is exactly solvable.

Modern computing power allows us to simulate highway and city traffic by looking at the microscopic behaviors of all cars. Perhaps the most well-known model of this kind is the cellular automaton based Biham-Middleton-Levine (BML) model of two-dimensional city traffic . Numerical simulations of this model strongly suggest that a first order phase transition from the full-speed phase to the completely jamming phase occurs as the car density of the system increases. Moreover, this phase transition is primarily the result of the exclusion volume effect . In spite of its simplicity, very little rigorous result is known for the BML model .

Different generations and variations of the BML model have been investigated in the literature . These generalizations focus on different aspects of the problem. They include the introduction of more realistic traffic rules , study of the effects of over-passes and faulty traffic lights , application of cellular automaton based traffic rules to single and multiple lane highway traffic  and investigation of higher dimensional traffic behaviors . In particular, Nagel and Schreckenberg  proposed a model of (one-dimensional) highway traffic. In addition to the regular acceleration and deceleration, they model the realistic behavior of car drivers by allowing them to apply their breaks in a stochastic manner. Later on, Fukui and Ishibashi (FI)  investigated a simple-minded deterministic model analogous to that of Nagel and Schreckenberg. Recently, Chowdhury and Schadschneider (CS) incorporated the one-dimensional highway traffic model of Nagel and Schreckenberg as well as the two-dimensional BML model together to study microscopic dynamics of city traffic .

It is not uncommon for local governments to set up tollbooths on highways. In fact, tolls can be used to control traffic flow and to increase revenue for local governments. Nevertheless, the presence of tollbooths will definitely slow down highway traffic. This is particularly true when tollbooths are set up on the highway rather than behind the entrance and exit ramps. Unfortunately, because of geographical and administrative considerations, tollbooths in many highways have to be built right on highways themselves. In fact, a few cellular automaton models of one-dimensional highway traffic flow with different kinds of blockages have been proposed and investigated . They consider the effects of overtaking sites, bottleneck and quenched noise. In contrast, this paper investigates how the presence of tollbooths affects the traffic flow in a single lane highway using a simple-minded deterministic discrete model based on the cellular automaton traffic model of Fukui and Ishibashi  as well as the so-called green wave model of Torok and Kertesz .

In their original model, Fukui and Ishibashi consider a one-dimensional array of $`N`$ sites with periodic boundary conditions. Each site may either be empty or have a single rightward moving car. They fix an integer $`V_{\mathrm{max}}`$ known as the maximum intrinsic car speed. Since their model does not consider the effect of car acceleration and deceleration, so at each timestep a car moves $`k`$ steps to the right where $`k`$ is the minimum of $`V_{\mathrm{max}}`$ and the number of consecutive empty sites to the left of the car. In addition, the motion of cars are updated in parallel . They define the average car speed by

$$V=\frac{1}{N\rho }\underset{i=1}{\overset{N\rho }{}}V_i,$$
(1)
where $`\rho `$ is the car density in the system (and hence $`N\rho `$ is the total number of cars in the system) and $`V_i`$ is the speed of the $`i`$th car. Note that the right hand side of Eq. (1) is averaged both over time and initial system configurations.

To coarse gain the highway system, we combine the FI model  with the so-called green wave traffic model . That is to say, we call a collection of consecutive sites all containing rightward moving cars a car cluster. We demand that cars in the same cluster to move altogether as a group with speeds equal to that of the leading car in the group except possibly when the car passes through a tollbooth. Besides, all updates are taken in parallel. Finally, we introduce the effects of tollbooths by selecting $`N_{\mathrm{booth}}`$ special sites on the system. Whenever a car reaches these special tollbooth sites, it has to stop immediately and to wait for $`t_{\mathrm{wait}}`$ timesteps before it is allowed to move again. This waiting car, therefore, may block the motion of the cars queuing behind it for $`t_{\mathrm{wait}}`$ timesteps. For simplicity, tollbooths are located uniformly on the system.

For example, if $`t_{\mathrm{wait}}=0`$, then a passing by car has to stop at the tollbooth site at once and then that car may move in the next timestep. Let us denote $`0`$ as a empty site, $`1`$ as a site occupied by a car, and an underline as a tollbooth site, then under our traffic rule, $`1100101\underset{¯}{1}00`$ will be transformed to $`0011010\underset{¯}{1}01`$ in the next timestep if $`V_{\mathrm{max}}=2`$, $`t_{\mathrm{wait}}=0`$.

Since we are interested only in the average car speed over all possible initial car configurations in the thermodynamic limit, therefore only the recurrent behavior of the system will affect the average car speed $`V`$. More precisely, $`V`$ depends only on the car density $`\rho `$, the maximum car speed $`V_{\mathrm{max}}`$, the tollbooth density $`d_{\mathrm{booth}}N_{\mathrm{booth}}/N`$ and the tollbooth stopping time $`t_{\mathrm{wait}}`$ of the system. Thus, $`t_{\mathrm{wait}}`$ and $`d_{\mathrm{booth}}`$ are the two controlling parameters in studying the effect of tollbooths.

In this paper, we are going to show that

###### Theorem 1

Case (a): If $`t_{\mathrm{wait}}=0`$, then the average car speed $`V`$ is given by

$$V=\{\begin{array}{cc}\left(d_{\mathrm{booth}}\frac{1}{d_{\mathrm{booth}}V_{\mathrm{max}}}\right)^1\hfill & \text{for }0<\rho \rho _a,\hfill \\ 1/2\rho \hfill & \text{for }\rho _a<\rho 1/2,\hfill \\ 1\hfill & \text{for }1/2<\rho <1,\hfill \\ 0\hfill & \text{for }\rho =1,\hfill \end{array}$$
(3)
where $`\rho _a=d_{\mathrm{booth}}1/d_{\mathrm{booth}}V_{\mathrm{max}}/2`$.

Case (b): If $`t_{\mathrm{wait}}>0`$, and $`1/d_{\mathrm{booth}}=1modV_{\mathrm{max}}`$ then the average car speed is given by

$$V=\{\begin{array}{cc}\left[d_{\mathrm{booth}}\left(t_{\mathrm{wait}}+\frac{1}{d_{\mathrm{booth}}V_{\mathrm{max}}}\right)\right]^1\hfill & \text{for }0<\rho \rho _b,\hfill \\ \left[\rho \left(t_{\mathrm{wait}}+1\right)\right]^1\hfill & \text{for }\rho _b<\rho <1,\hfill \\ 0\hfill & \text{for }\rho =1,\hfill \end{array}$$
(4)
where $`\rho _b=d_{\mathrm{booth}}(t_{\mathrm{wait}}+1/d_{\mathrm{booth}}V_{\mathrm{max}})/(t_{\mathrm{wait}}+1)`$.

Case (c): For the remaining possibility, namely that $`t_{\mathrm{wait}}>0`$ and $`1/d_{\mathrm{booth}}1modV_{\mathrm{max}}`$,, the average car speed is given by

$$V=\{\begin{array}{cc}\left[d_{\mathrm{booth}}\left(t_{\mathrm{wait}}+\frac{1}{d_{\mathrm{booth}}V_{\mathrm{max}}}\right)\right]^1\hfill & \text{for }0<\rho \rho _{c1},\hfill \\ \left[\rho \left(t_{\mathrm{wait}}+2\right)\right]^1\hfill & \text{for }\rho _{c1}<\rho \rho _{c2},\hfill \\ \left[d_{\mathrm{booth}}\left(t_{\mathrm{wait}}+1+\frac{1}{d_{\mathrm{booth}}V_{\mathrm{max}}}\right)\right]^1\hfill & \text{for }\rho _{c2}<\rho \rho _{c3},\hfill \\ \left[\rho \left(t_{\mathrm{wait}}+1\right)\right]^1\hfill & \text{for }\rho _{c3}<\rho <1,\hfill \\ 0\hfill & \text{for }\rho =1,\hfill \end{array}$$
(5)
where $`\rho _{c1}=d_{\mathrm{booth}}(t_{\mathrm{wait}}+1/d_{\mathrm{booth}}V_{\mathrm{max}})/(t_{\mathrm{wait}}+2)`$, $`\rho _{c2}=d_{\mathrm{booth}}(t_{\mathrm{wait}}+1+1/d_{\mathrm{booth}}V_{\mathrm{max}})/(t_{\mathrm{wait}}+2)`$, $`\rho _{c3}=d_{\mathrm{booth}}(t_{\mathrm{wait}}+1+1/d_{\mathrm{booth}}V_{\mathrm{max}})/(t_{\mathrm{wait}}+1)`$.

Before going on to prove this theorem, we remark that in all the three cases above, first order phase transitions in $`V`$ occur only at $`\rho =1`$. All other transition points are second order in nature. (See Fig. 1a–c for typical shapes of the $`\rho `$ versus $`V`$ curves in these three cases.)

Now, we begin to prove this theorem by first considering the following Lemmas and Propositions.

###### Proposition 1

The speed of cars in a recurrent state is never limited by another car provided that the car density of the system $`\rho `$ is less than or equals to $`\rho _a`$, $`\rho _b`$ and $`\rho _{c1}`$ for cases (a), (b) and (c) in Theorem 1 respectively. In other words, cars in systems with the above parameters move in maximum possible speed except possibly when they enter tollbooth sites.

Proof: According to our traffic rules, cars pass through a tollbooth site at a rate of at most one car every $`(t_{\mathrm{wait}}+1)`$ timesteps. Moreover, a car cannot overtake one another and a car cluster can only break up at a tollbooth site. Thus, if there is a car in the recurrent state whose speed is limited by another car in front (towards the downstream direction), then all cars in the recurrent state must be blocked by the car in front therefore preventing them from moving in maximum speed $`V_{\mathrm{max}}`$ some time under the repeated application of our traffic rules. More importantly, the number of timesteps required to travel from one car to the next in a recurrent state under the action of our traffic rule cannot exceed $`t_{\mathrm{wait}}+1`$. And, the blocking car must in turn being blocked in at least one of the previous $`(t_{\mathrm{wait}}+1)`$ timesteps. Inductively, this can happen only when the average distance between cars in the configuration is greater than the minimum possible distance of a full speed recurrent configuration. But in a full speed recurrent configuration, a car can only be found at a tollbooth site, or at a site $`kV_{\mathrm{max}}`$ steps ahead of a tollbooth site where $`k=1,2,3,\mathrm{},1/d_{\mathrm{booth}}V_{\mathrm{max}}`$ and a car takes $`(t_{\mathrm{wait}}+1/d_{\mathrm{booth}}V_{\mathrm{max}})`$ timesteps to move through two consecutive tollbooths. Thus, there are $`(t_{\mathrm{wait}}+1/d_{\mathrm{booth}}V_{\mathrm{max}})`$ distinct possible “locations” for the cars in a recurrent state between two successive tollbooths if cars waiting at the same tollbooth for different times are regarded as distinct locations. And in order that cars do not block one another, an average of at most $`0.5`$, $`\mathrm{min}(0.5,1/(t_{\mathrm{wait}}+1))=1/(t_{\mathrm{wait}}+1)`$ and $`\mathrm{min}(0.5,1/(t_{\mathrm{wait}}+2))=1/(t_{\mathrm{wait}}+2)`$ of these “locations” are occupied in cases (a), (b) and (c) respectively. (Note that the difference between case (b) and (c) is due to the fact that in case (b) an unobstructed car can move into an occupied tollbooth site at the same time when that occupying car moves out of that tollbooth site for $`1/d_{\mathrm{booth}}=1modV_{\mathrm{max}}`$. In contrast, this can never happen in case (c).) Hence, in order that some car in the recurrent state is being blocked by another car, $`\rho `$ must be greater than $`\rho _a`$, $`\rho _b`$ and $`\rho _{c1}`$ in case (a), (b) and (c) respectively. $`\mathrm{}`$

###### Lemma 1

Provided that $`t_{\mathrm{wait}}>0`$ and $`V_{\mathrm{max}}>1`$, a car in a recurrent state can only be blocked by a car cluster passing through a tollbooth.

Proof: Suppose the contrary, we can find car A blocking car B in a recurrent state such that car A does not belong to any car cluster passing through a tollbooth. If we denote the action of applying our traffic rule once by $`𝐓`$, then the inverse map $`𝐓^1`$ is well-defined on the set of all recurrent states. Clearly, under the action of $`𝐓^1`$, car A must move backward by at least one site. Moreover, car A must be blocked by another car at least once in the previous $`(t_{\mathrm{wait}}+1)`$ timesteps. Inductively, by considering the repeated application of $`𝐓^1`$, we end up with a car C located at a tollbooth site whose motion in the next timestep is blocked by a car D located two sites in front. But this is impossible as the inverse image of this configuration under $`𝐓^2`$ is an empty set, contradicting the assumption that the state is recurrent. $`\mathrm{}`$

###### Proposition 2

In case (c), whenever the car density lies between $`\rho _{c1}`$ and $`\rho _{c2}`$, a car in a recurrent state can only be blocked by another car right at a tollbooth site when the car occupying that tollbooth site is moving away. And whenever the car density lies between $`\rho _{c2}`$ and $`\rho _{c3}`$, a car in a recurrent state can only be blocked by another car right at a tollbooth site when the car occupying that tollbooth site is going to move away in the next timestep.

Proof: Proposition 1 implies that cars begin to block one another when $`\rho >\rho _{c1}`$ in case (c). From Lemma 1 and the fact that $`1/d_{\mathrm{booth}}1modV_{\mathrm{max}}`$, we know that it is possible for a car A in a recurrent state to be blocked by a car located at a tollbooth site. If this event happens, car A has to take $`(t_{\mathrm{wait}}+1+1/d_{\mathrm{booth}}V_{\mathrm{max}})`$ timesteps to move through that two successive tollbooth sites. In other words, car A takes one timestep more than the minimum possible value in order to move through the two successive tollbooths. Note that if $`\rho \rho _{c3}`$, no recurrent state can contain a car cluster C making up of more than two cars. The reason is simple: for otherwise, there exits, at any instance, at least one interval J between two tollbooths containing no more than $`(t_{\mathrm{wait}}+1/d_{\mathrm{wait}}V_{\mathrm{max}})/(t_{\mathrm{wait}}+2)`$ cars. But then the car outflow rate from this interval J is strictly less than one car per $`(t_{\mathrm{wait}}+1)`$ timesteps while from Lemma 1 the car outflow rate from the cluster C equals one car every $`(t_{\mathrm{wait}}+1)`$ timesteps. Thus eventually there is not enough car supply to maintain the car cluster C and new car clusters with more than two cars cannot be found elsewhere in the system due to the restrictions of both the can inflow and outflow rates in a tollbooth site. This contradicts our assumption that the configuration is recurrent.

Using the same trick as in the above argument that no car cluster with more than two cars can be formed in a recurrent configuration for $`\rho \rho _{c3}`$, it is easy to show that when $`\rho (\rho _{c1},\rho _{c2}]`$, the recurrent state consists of intervals of freely moving cars with density $`\rho _{c1}`$ as well as intervals of cars with density $`\rho _{c2}`$ that can be blocked by a ready-to-move car in a tollbooth site. (All car densities mentioned here are averaged over $`(t_{\mathrm{wait}}+2)`$ timesteps.)

Similarly, when $`\rho _{c2}<\rho \rho _{c3}`$, a recurrent state is made up of intervals of cars with density $`\rho _{c2}`$ that can be blocked by a ready-to-move car in a tollbooth site as well as intervals of cars with density $`\rho _{c3}`$ that can be blocked by a car in a tollbooth site that will move in the next timestep. $`\mathrm{}`$

###### Lemma 2

Let $`t_{\mathrm{wait}}=0`$ and $`V_{\mathrm{max}}>1`$. Then if $`\rho 1/2`$, all car clusters in a recurrent state consist of only one car. And if $`\rho 1/2`$, the recurrent state contains no consecutive empty site.

Proof: Since $`t_{\mathrm{wait}}=0`$, all cars can move at least one step to the right in every timestep. Besides, it is not possible for two car clusters to merge. Since we are using green wave traffic rule and $`V_{\mathrm{max}}>1`$, a car cluster may break up into two only at a tollbooth site. Furthermore, consecutive empty sites will “move” to the left while car clusters move to the right. Thus in $`\text{O}(N)`$ timesteps, the system will evolve to a state with maximum possible number of car clusters. Hence the lemma is proved. $`\mathrm{}`$

###### Proposition 3

Let $`t_{\mathrm{wait}}=0`$, $`V_{\mathrm{max}}>1`$ and $`\rho _a\rho <1/2`$, then exactly one car passes through a tollbooth every two timesteps.

Proof: From Lemma 2 and our traffic rule, we know that at most one car can pass through a tollbooth every two successive timesteps. And from Proposition 1, we know that at least one car is blocked in a recurrent state. Using a similar argument as in the proof of Proposition 2, one can always find a car whose motion is blocked by another about-to-go car locating at a tollbooth site in front. Hence, every car takes exactly two timesteps to pass through such a tollbooth. $`\mathrm{}`$

###### Proposition 4

The average car speed formulae in case (b) and (c) are valid for $`\rho _b<\rho <1`$ and $`\rho _{c3}<\rho <1`$, respectively.

Proof: From Lemma 1 and the assumption that $`\rho _b<\rho <1`$ and $`\rho _{c3}<\rho <1`$ in case (b) and (c) respectively, we can always find a car cluster in a recurrent state making up of at least three cars lining up in front of a tollbooth. Moreover, such a car cluster can never dissolve (that is, the cluster never disintegrates to clusters of single cars) under the repeated action of our traffic rule. Since the leading car in this car cluster passes through the tollbooth at a rate of once every $`(t_{\mathrm{wait}}+1)`$ timesteps, we conclude that the average car speed $`V`$ equals $`1/\rho (t_{\mathrm{wait}}+1)`$. $`\mathrm{}`$

Proof of Theorem 1: Since $`V`$ is clearly equal to 0 when $`\rho =1`$, so it remains for us to consider the case when $`\rho <1`$.

Let us first consider case (a). If $`V_{\mathrm{max}}=1`$, then our traffic rules reduce to moving each car forward one site at a time provided that $`\rho <1`$. Hence, Eq. (3) is trivially true in this case. So we only need to consider the case when $`V_{\mathrm{max}}>1`$. And from Proposition 1, we know that for case (a), when $`\rho \rho _a`$, cars will eventually moves as if there were no other cars in the system. Thus a car in a recurrent state takes $`(t_{\mathrm{wait}}+1/d_{\mathrm{booth}}V_{\mathrm{max}})=1/d_{\mathrm{booth}}V_{\mathrm{max}}`$ timesteps to travel through $`1/d_{\mathrm{booth}}`$ sites. Hence, $`V=\left(d_{\mathrm{booth}}1/d_{\mathrm{booth}}V_{\mathrm{max}}\right)^1`$. And when $`\rho _a<\rho 1/2`$, $`V=1/2\rho `$ is an immediate consequence of Proposition 3. Finally, when $`1/2<\rho <1`$, Lemma 2 tells us that two adjacent car clusters are separated by exactly one empty site. Hence, $`V=1`$ and Eq. (3) holds for $`t_{\mathrm{wait}}=0`$.

Now we are going prove Eq. (4) in case (b). Clearly the validity of Eq. (4) for $`\rho \rho _b`$ follows directly from Proposition 1. And finally $`V`$ for $`\rho _b<\rho <1`$ has just been proven in Proposition 4.

Lastly, the validity of Eq. (5) in case (c) when $`\rho \rho _{c1}`$ or $`\rho >\rho _{c3}`$ follows directly from Propositions 1 and 4. And if $`\rho _{c1}<\rho \rho _{c2}`$, Proposition 2 tells us that car density (averaged over $`(t_{\mathrm{wait}}+2)`$ timesteps) between any two successive tollbooths in a recurrent state is either equal to $`\rho _{c1}`$ or $`\rho _{c2}`$. Thus, with a steady and continual supply of cars from behind, the car density averaged over $`(t_{\mathrm{wait}}+2)`$ timesteps in every interval between two successive tollbooths is a constant. Since a car is released once every $`(t_{\mathrm{wait}}+2)`$ timesteps from the tollbooth in a an interval with car density $`\rho _{c2}`$, hence $`V=1/\rho (t_{\mathrm{wait}}+2)`$ in this car density range. Finally, when $`\rho _{c2}<\rho <\rho _{c3}`$, Proposition 2 implies that cars pass through every tollbooth once every $`(t_{\mathrm{wait}}+2)`$ timesteps. Hence, Eq. (5) holds in this density range as well. $`\mathrm{}`$

In summary, we have investigated the behavior of a single lane deterministic highway traffic model in the presence of tollbooths. Our models are exactly solvable and the average car speed consists of a high speed, a partially jamming and a trivial completely jamming phases. The transition from the partially jamming phase to the completely jamming phase is first order in nature while all other transitions are second order.

Most importantly, our model suggests that the average car speed at high car density depends only on the car density and is independent of the detail arrangement of tollbooths in single lane traffic. While the regular placement of tollbooths and the deterministic traffic rules give rise to the unrealistically flat $`V`$ when car density lies between $`\rho _{c2}`$ and $`\rho _{c3}`$, the general observation that for $`t_{\mathrm{wait}}>0`$, the average car speed $`V`$ is approximately inversely proportional to $`t_{\mathrm{wait}}`$ is robust. Let us compare our results with the two-dimensional CS city traffic model. In the CS model, the time duration of red or green lights plays an analogous role of the tollbooths. As shown in their density versus flux per street curve in Fig 4 of Ref. , a linear region corresponding to cars moving with almost full speed is observed when the car density $`\rho `$ is low. More interestingly, a plateau region corresponding to the $`V1/\rho `$ is observed when the time duration of red or green lights is large and $`\rho `$ is around 0.1 to 0.4. And $`V`$ starts to decrease at higher values of $`\rho `$ probably due to the effects of two-dimensional car blocking. Comparing the observations in the CS and our models, we believe that the $`1/\rho `$ behavior in high car density highway traffic is robust.

###### Acknowledgements.

H.F.C. and H.X. are supported in part by the Hong Kong SAR Government RGC Grant HKU 7098/00P and H.F.C. is also supported in part by the Outstanding Young Researcher Award of the University of Hong Kong.