text:
# Dynamics of polymer chain collapse into compact states

## I Introduction

Several decades of protein folding simulation have produced a substantial body of knowledge about processes governing molecular collapse from a random denatured state to a well-ordered and uniquely defined ground state Brooks III et al. (1988); Dill et al. (1995); Shakhnovich (1997); Brooks III (1998); Pande et al. (2000). Simulations addressing detailed representations of relatively large molecules over extended time periods lie close to, or beyond, the capabilities of currently available computers; consequently, a variety of simplified models designed to capture particular aspects of the molecular behavior have been formulated. These models eliminate much, or even most of the molecular detail, thereby achieving a major reduction in the amount of computation required for evaluating the interactions that govern the behavior. Further reduction in effort is gained by replacing Newtonian dynamics by one of several forms of Monte Carlo sampling procedure, and the configuration space available to the molecule is generally reduced substantially by discretizing the problem and confining the molecule to the sites of a regular lattice Socci and Onuchic (1994); Dill et al. (1995); Dinner et al. (1996). After imposition of these and other approximations the connection between the models and the original problem is at best tenuous, so that it is not always obvious if and how conclusion arising from such models relate to real proteins.

Although a full dynamical simulation of a relatively detailed representation of even a small protein (or part of a protein) remains a major computational undertaking Duan and Kollman (1998), there is no reason why simplified models cannot be studied in the continuum with molecular dynamics (MD) methods. The purpose of the present paper is to extend an earlier study of this kind Rapaport (2002), that dealt with helix formation, to address the formation of a series of more complex native structures; the expectation is that through such simulations information will emerge about how structural variation influences the folding process.

The native structure chosen for study is the cube. The molecule itself is represented as a linear chain, free to move in the continuum (subject to specified internal degrees of freedom), and the interactions are chosen in such a way that, in the native low-energy state, the chain is folded so that its sites form a perfectly ordered cubic array. There are numerous ways that a linear path can wend its way through all the sites of a cubic lattice of a given size; from this collection, four distinct but conveniently described paths have been selected for study. Cubic conformations are not known to occur in real proteins, but this is not an issue since the purpose of the study is to examine the effect of the variation of ground-state configuration on folding behavior. As will become apparent subsequently, the difference between these states directly impacts both the likelihood of entanglement (misfolding) and the rate of convergence to the ordered state (correct folding). This difference is at least partly attributable to the dissimilar contact patterns: adjacent sites in the native state involve different proportions of nearby and more distant (or local and nonlocal) sites when measured in terms of backbone separation.

In dealing with highly detailed protein models, there is a problem establishing that the free-energy minimum of the overall potential function really corresponds to the known native conformation; given the present state of the art this is practically impossible, owing to the way potentials are developed and the limited spatial resolution of experimental structure determination. The alternative approach, based on a highly simplified description, preserves the underlying tenet of protein structure theory, namely that the primary sequence determines structure, but is sufficiently economical to allow complete folding pathways to be studied. Moreover, the computational efficiency allows the dynamics of not just a single folding process to be followed, but an entire ensemble of systems can be examined, leading to a representative sample of the kinds of behavior that can arise; the importance of extensive sampling cannot be overstated, since with a very small number of samples it is impossible to determine what is in fact “typical” behavior.

One reason for choosing a cubic form for the native state is that a considerable body of work exists on chains that fold into this overall shape (see Dinner et al. (1996) and references therein). The difference is that all such work involves lattices, primarily using various Monte Carlo techniques; calculations of this type are carried out at a constant nominal temperature (the temperature is a parameter of the Monte Carlo procedure used in deciding whether to accept randomly generated trial configurations, and it determines the ability to cross potential-energy barriers). In lattice studies, the potential energy is a sum over contributions from chain sites that are in contact because they occupy adjacent lattice sites; the approach to choosing the potential and deciding which site pairs should be encouraged to form contacts depends on the model, and typically the interactions are chosen in such a way that the native state has a substantially lower potential energy than any other state, even those corresponding to different cube packings (in some of these studies, interactions are selected so that even in the ground state a fraction of adjacent site pairs would prefer not to be neighbors, thereby allowing a certain amount of frustration to be incorporated into the design). The folding behavior has been found to depend on the gap between this global energy minimum and the other local minima corresponding to compact but misfolded configurations. Much of the lattice work has dealt with cubes having three sites per edge, corresponding to relatively short chains of length 27. When chains of length 125 were studied by similar means Dinner et al. (1996) it was noted that results based on shorter chains were subject to finite-size effects; in particular, the energy gap was a necessary but not sufficient condition for folding, and the presence of a set of core sites, nonexistent in 27-mers, played a key role in the folding.

While it is possible to debate the merits of Monte Carlo relative to MD, there is at least one difference that is particularly important for folding studies. In lattice-based Monte Carlo, the spatial discretization imposes restrictions on the permissible internal rearrangements: configuration changes typically involve the displacement of a single chain site, or a very small number of adjacent sites, as in crankshaft-type motion, and all displacements are based on jumps between lattice sites. There is no provision for collective movement of more extended subunits, a mode of reorganization likely to dominate once a chain has reached an even moderately compact state, and clearly such restrictions will have an impact on the chain “dynamics”. The MD approach is free from any limitations of this kind and, at least in principle, aims at a realistic representation of the dynamics.

The goal of the present paper is to extend the approach used previously Rapaport (2002) for modeling helix collapse to the study of cubic configurations, with emphasis on how the folding process is affected by the way the chain is arranged inside the cube. The main focus of the earlier paper was on the ability of a chain having a known native state to actually fold into that structure, within the course of a simulation of reasonable duration; the obvious extension of the work is to consider other configurations with different structural features. A helix involves contacts between sites relatively closely spaced along the chain, whereas cubic structures involve varying mixtures of local and nonlocal neighbors, a distinction that is likely to affect the folding process. The following sections outline the techniques, insofar as they differ from the earlier work, and then proceed to a discussion of the results and their implications.

## II Methodology

The model follows the approach described previously Rapaport (2002) in which the chain is represented as a series of sites (or masses) linked by bonds of constant length, and the angles between successive bonds are assigned fixed values consistent with the desired native ground state; the only internal degrees of freedom are the dihedral angles that describe the relative orientations of next-neighbor bond pairs projected in a plane perpendicular to the bond between them. The interactions included in the model involve a soft-sphere repulsion between sites (the short-range, repulsive part of the Lennard-Jones potential) that is responsible for the excluded volume of the chain, and a torsional potential acting along each bond (except for the first and last bonds) whose minimum corresponds to the dihedral angle value in the native state. Aside from the different bond and dihedral angles as described here, all sites are identical. While both of these interactions appear in more realistic protein potentials, the present model does not include any of the site-specific pair interactions that are principal component of potentials used in more detailed protein studies; since there are no direct interactions between sites, any contact preferences exhibited by individual sites are entirely due to the torsional interactions that act along the backbone. Additionally, since the presence of a solvent would slow the computations substantially, the chain resides in a vacuum (a characteristic of some more detailed simulations as well); this also avoids the need to decide on the level of detail used to describe solvent effects. At a qualitative level, omissions such as these should not prevent simplified models from providing insight into the underlying dynamics of folding.

The MD approach used for the simulation involves recursive techniques for dealing with the internal degrees of freedom, while the integration of the equations of motion (for translational, rotational and internal motion) is based on the leapfrog algorithm; the technical aspects of computations of this kind have been addressed at length elsewhere Rapaport (2002, 2003) so that the details need not be repeated here. The chain is initially heated to comparatively high temperature to produce a random configuration, and is then cooled very slowly by reducing the kinetic energy by a constant factor at regular time intervals. The first stage of the cooling process is designed to extract energy from the chain and leave it trapped in a free-energy well – preferably, though not always, one from which the native state is accessible – from which escape is unlikely; the second stage is intended to freeze the configuration by eliminating, or at least weakening, the soft modes that correspond to the low-frequency collective motion present even at low temperature, so that measurements of the properties of the configuration eventually reached are relatively free from the effects of thermal fluctuation.

As in the helix studies Rapaport (2002), the cooling rate is chosen empirically, and is dependent on chain length. In order to allow comparison of the different native-state chain arrangements, the same cooling rate is used in all cases. The rate value is chosen to differentiate clearly between the varying levels of ability to find the correct folding pathway to the native states; cooling too fast will allow none of the chains to fold properly, while excessively slow cooling will result in a situation where there are practically no unsuccessfully folded chains.

As an alternative (not explored here) to gradual cooling, the simulation could be carried out at constant temperature, using appropriate MD techniques, with a temperature value chosen to achieve a reasonably rapid folding rate while avoiding the excessive entanglement responsible for misfolding. Such results would produce a long-tailed (or open-ended) distribution of folding times, so that ensuring adequate configuration sampling is likely to demand much longer computation times than the approach based on cooling. The information produced concerning the comparative folding abilities of different native configurations and the behavior observed along the folding pathways should be similar, although timescales will differ due to the way temperature is involved.

The initial angular velocities associated with the dihedral angles are randomly assigned; changes to the seed used for initializing the random number generator lead to entirely different folding scenarios. Since it is neither possible, nor particularly meaningful, to describe at length the detailed histories of each of the large number of runs carried out, time-dependent ensemble averages over the collection of histories are evaluated. In some of the analyses described below the results are divided into two groups, runs that produce successful folding and runs that become trapped in a misfolded state; other kinds of analysis apply to either the entire set of runs or just the successful folders. While the results could be divided according to other criteria, that based on successfully completed folding is the most directly relevant.

Chains of length 64 and 125 have been studied; these chains can be packed into cubic states with four and five sites per edge, respectively. In terms of the reduced unit of length, chosen here (as usual) to be the parameter $`\sigma `$ of the Lennard-Jones interaction (which, for the soft-sphere case, is slightly less than the sphere’s effective diameter) the length of the fixed links between neighboring chain sites is 1.3; this bond length is sufficiently short to prevent self-intersection, while allowing physically larger molecules than if adjacent spheres were actually touching. A chain of $`N`$ sites has $`N1`$ links and, since torsional motion is not associated with the two end links, $`N3`$ internal degrees of freedom. Bond angles are fixed at either $`\pi /2`$ (right angle) or $`0`$ (parallel), depending on the location in the native-state cube; to ensure that all bonds are treated equivalently, torsional motion is associated with bonds even when successive bonds are in the same direction. The torsional potential is the same sinusoidal function used in the helix studies, with the energy minimum located at the native-state dihedral angle for each individual bond.

The runs are of length $`10^6`$ and $`1.5\times 10^6`$ time steps for the $`N=64`$ and $`N=125`$ chains respectively; the size of the time step is 0.004 in reduced units. Cooling is accomplished by scaling the instantaneous kinetic energy every 4000 steps by factors of 0.97 or 0.98 for the shorter and longer chains, a larger factor implying slower cooling. This results in a kinetic energy fall from an initial value of approximately 4 (per degree of freedom, in reduced units defined in terms of the soft-sphere potential) to a final value of $`2\times 10^3`$ (the sum of potential and kinetic energy first becomes negative about a tenth of the way through the run); by comparison, the potential energy per dihedral angle is $`5`$ in the ground state. A total of 400 runs are carried out for each type of native state and chain length; measurements and configurational snapshots are recorded at suitable time intervals.

The four kinds of chain organization used for the native states are shown in Fig. 1; here, for clarity, the chains appear as continuous tubes, whereas in actual fact each chain consists of a series of closely-spaced spheres joined by bonds of fixed length (an example of this representation is shown subsequently). In the orientation shown, the cubes are filled one complete horizontal plane at a time; the configurations shown are for $`N=125`$, the organization for $`N=64`$ chains is similar, after allowing for the smaller cube size. The differences between the configurations relate to the manner of filling the planes and the relationship between adjacent planes. Two of the cases, labeled A and B, are based on zigzag (or antiparallel) patterns, the other two, C and D, are spirals. The other difference is the relation between the fill patterns of successive planes: the two zigzag configurations consist of rotations (A) and reversals of the fill pattern (B); the two spiral patterns consist of spirals that rotate in alternate directions (C) and spirals that continue to rotate in the same direction (D). In all cases, the chain terminal sites lie on the outer faces of the cube. There are numerous other possible patterns (some more readily characterized than others), in particular, patterns that are not restricted to filling a single plane at a time (such as the three-dimensional Hilbert curve), but the present assortment is already sufficient for displaying a wide range of behavior. The filling sequences create a distinction between “secondary” (each plane) and “tertiary” (the entire cube) structural features which some alternative pathways through the cube might not exhibit. Fig. 2 shows an early, essentially random state from one of the runs after the chain has been heated to erase memory of the initial state; here a ball-and-stick representation is used to show the actual chain design.

## III Results

The analysis focuses on a selection of configurational averages and distributions that characterize the time-dependent behavior of the chains as folding progresses. Some quantities will turn out to be sensitive to the organization of the native state, namely the way the cube is packed, others less so. Although averaging merges the details of individual folding histories, it does facilitate the extraction of key aspects of the behavior that are contained in what is sometimes widely scattered data; the spread of the distributions, where shown, provides an indication of the inherent variability of the behavior.

### III.1 Folding success rate

Ultimately, the only truly reliable test of whether a configuration has folded successfully is based on visual examination, but while this technique has been employed, it does not offer a reasonable approach for describing large sets of configurations. Furthermore, owing to residual thermal fluctuations, the separation between sites that are deemed to be within contact range will usually exceed the value applicable to a perfectly formed cube, and this flexibility must be accommodated by the analysis. A corresponding difficulty does not arise in lattice studies, where the existence of a contact requires the sites to be neighbors on the lattice in which the chain is embedded.

The automated decision scheme introduced here to classify the chain states employs two defining parameters. The first is a factor $`\varphi `$ multiplying the bond length that establishes the maximum contact range. Visual inspection suggests $`\varphi =1.3`$; this excludes essentially all false positives but can underestimate the success rate by omitting some slightly misaligned but otherwise correct configurations – small thermal vibrations (generating concertina-like modes that produce slightly nonparallel layers, or minor deviations from planarity within the layers) are permitted energetically since the cost is comparatively low. Raising the value to $`\varphi =1.4`$ improves the situation by capturing most of these configurations, but allows a small number of false positives. The second parameter is the minimal fraction of site pairs that must satisfy the distance criterion for the chain to be regarded as successfully folded. Here a value of 0.8 is used, and this fraction must appear in at least one configuration during the final 120 time units (amounting to a few percent) of the run; in practice, once this value is encountered, a high contact fraction is normally maintained throughout the remainder of the run.

The fractions of chains of each type that fold successfully over the entire set of runs are summarized in Table 1. The results are ranked in descending order; this order has also been used to assign labels to the native-state configurations. The success rates for both values of $`\varphi `$ are shown; the difference is either zero or quite small. Thus, for practical purposes, the results are insensitive to the choice of $`\varphi `$ in this range, and a value of $`\varphi =1.3`$ will be used when required in subsequent analysis.

For each configuration type the success rate is lower for the longer chains, but for both cube sizes there is a significant, systematic variation that depends on the organization of the native state. The most successful folder is type A, which consists (see Fig. 1) of successive zigzag layers in a crossed alignment. Relatively similar success rates are achieved by types B and C, which are the antiparallel (or reversed-direction) zigzag layers and the reversed-direction spirals; the common feature of both these configurations is that successive layers are mirror images of one another (leading to prominent features in the contact maps discussed later). Finally, the least successful folders are those of type D, in which the spirals, both increasing and decreasing, have all their turns in the same direction. The fraction of successful folders as a whole depends on the cooling rate and, as indicated earlier, the rates were chosen to allow differentiation between the different chain types; the effect of alternative cooling rates is demonstrated in Rapaport (2002) for the case of helix folding, and cubes respond in a similar manner.

Fig. 3 shows examples of misfolded states, one for each chain type; much of the correct native structure is present, but a very small number (as small as one) of incorrect twists leads to the wrong overall organization. Visual examination of numerous final states leads to the conclusion that while all correctly folded chains are alike, each incorrectly folded chain is incorrect after its own fashion. One could attempt to classify incorrect folds: in the helix study, most misfolds were due to a single incorrect twist, while here, an example of a clearly identifiable misfold involving just a single twist is the case in which an outer layer is rotated out of alignment and jammed against one of the other faces of the cube. However, it is not apparent that any useful purpose would be served by a systematic search among the misfolded states for common structural motifs. Residual thermal vibration was the reason given earlier for failing to identify states as correctly folded; an example of a folded configuration representing a borderline case is shown Fig. 4.

### III.2 Energy and radius of gyration

The most familiar global properties of flexible chain molecules are the internal energy and the radius of gyration (another quantity familiar from chain-configuration studies, the mean-square end-to-end separation, does not convey sufficient information to be useful for analyzing folding). The distance of each of these quantities from the known native-state value provides an averaged measure of the deviation of the current chain state from the ordered configuration, without any attempt at characterizing the precise form of this deviation. These quantities, unlike some more specific measures described later, also correspond to properties that, at least in principle, are measurable (by calorimetry and light scattering) in the laboratory.

The graphs of normalized internal energy (the only interactions incorporated in the model, apart from excluded volume which makes a very small contribution, are the torsional terms) are shown in Fig. 5 for $`N=125`$; since the energies are negative the magnitudes are shown. The overall average for each of the four sets of configurations is plotted as a function of time, with separate curves used to distinguish the averages for successful and unsuccessful folders (using the definition of success given above). It is clear from the graph that there is very little to separate the different sets of data, implying that energy is not a useful distinguishing factor for these studies. Examination of the full energy distributions (not shown) reveals the spread in values to be very narrow. As will become apparent subsequently, the internal energy approaches the limiting value of the folded state well before other quantities achieve convergence.

An overall measure of the compactness of the configuration is provided by the radius of gyration; it is defined as the root-mean-square distance of sites (assuming all to have the same mass) from their mutual center of mass (if the shape deviates significantly from spherical – or cubical – symmetry, then it might be necessary to consider the individual moments of the mass distribution, but this is not the case here). The mean radius of gyration is shown in Fig. 6, again with separate curves for successful and unsuccessful folders of each configuration type. The results are normalized relative to the correctly folded collapsed state, which for a bond length of 1.3 has the value 10.14 for a packed cube of size five (6.34 for size four). The results exhibit some degree of spread over much of the run duration, reflecting different collapse rates, but the final results are very closely spaced.

A common feature of the results in Fig. 6, for all configuration types but more pronounced in some cases, is that the unsuccessful folders appear to collapse somewhat more rapidly during the initial stage of folding, only to be overtaken later on by the chains that do manage to fold correctly. There is very little difference in the final mean values for the different configurations and folding outcomes. The actual distributions of values (not shown) are fairly narrow, but do reveal additional details; for example, the distribution at the end of the runs for type C configurations consists of a pair of narrow partially overlapping peaks, although the combined peak width is no greater than for type D with only a single peak. Since it follows from Fig. 6 that the range of actual (unnormalized) mean values is less than the site diameter it is apparent that, as with the internal energy, the radius of gyration is of little help in differentiating between sets of correctly and incorrectly folded states (the more extreme cases of misfolding make only a small contribution to the mean). Though the radius of gyration converges less rapidly than the energy, it is already close to its limiting value at a time when other chain properties (see below) are still undergoing substantial change. This initial collapse to a relatively compact state, followed by a (sometimes successful) search for the route to the native folded state, is also a feature of lattice Monte Carlo studies Dinner et al. (1996); however, given the continuity of configuration space in the MD approach, it is likely that MD will be more effective in rearranging a chain once it has reached a semi-compact state.

### III.3 Contact formation

The study of properties based on the local environments of the chain sites provides an alternative to the global measures of energy and conformation considered above. An order parameter $`S`$ that measures the proximity of a configuration to its native state can be defined in terms of the fraction of site pairs forming contacts in the native state that lie within contact range – with the separation factor $`\varphi `$ introduced earlier used to determine which pairs qualify (for a chain whose native state is a cube with $`n`$ sites per edge, the maximum number of contacts is the number of lattice bonds $`n(n1)(3n2)`$ minus the number of chain bonds $`n^31`$). The time dependence of $`S`$ is shown in Fig. 7. Here, and subsequently, surface plots are used to display these time-varying distributions in compact fashion. The fraction of successful folders corresponds to the area under the peak in the region close to unity near the right corner of the grid; unsuccessful folders contribute to lower values. Clearly, the behavior of $`S`$ depends not only on chain length but also on the cube packing. In some cases, the $`S`$ distribution becomes bimodal in the course of the folding process, with the two peaks corresponding to chains that do and do not manage to fold successfully ($`S0.8`$ is the minimal value exhibited by a folded chain, in which a fraction of the pairs that should be classified as being in contact is temporarily out of range due to thermal fluctuations); in other instances, a single broad peak encompasses both folded and misfolded configurations, and its position depends on the configuration type.

An alternative approach to describing the degree of contact formation that avoids the need for the factor $`\varphi `$ is based on measuring the separation distribution of site pairs that should lie within contact range. The results are shown in Fig. 8; in those cases where multiple peaks occur, it is the narrow peak located at the lowest value that corresponds to correctly positioned neighbors, with peaks at larger separations showing the distance distribution of those site pairs that are unable to position themselves correctly; the graphs include all chains, both folders and non-folders, with both kinds contributing to the various peaks, although in different proportions.

Contact formation between particular site pairs can be selectively studied by classifying the sites into different categories and evaluating the average behavior for each category. One such classification, used previously in lattice studies Dinner et al. (1996), distinguishes sites belonging to the exterior surface faces of the cube from those of the interior core. Although the design of the native states of the chains considered here does not assign any preferred role to such a core, it does not preclude analysis of the results from such a perspective. Chains that fill cubes of size four have a mere eight sites in the core (12.5% of the sites), but chains filling cubes of size five have a more substantial core of 27 sites (21.6%). Fig. 9 shows the contact fractions for chains that are able to fold successfully (the two sets of contact distributions – for core and surface sites – are normalized separately and then combined in each plot); with the longer chains, it is apparent that for configurations C and D contacts form, on average, at different rates and the core sites tend to make contact with their correct neighbors earlier (separate plots of the distributions establishes this order). Such behavior is indicative of structure (in an averaged sense) nucleating in the interior and propagating outwards, but it is not apparent for all chain types and lengths.

An alternative to making the distinction between core and surface sites is to classify the individual contacts according to their role in the structure. All the chain configurations considered in the present work share the common design characteristic of a native structure built from a series of planar layers. This suggests dividing the contacts into two classes, those between sites within the same layer (each layer can be regarded as a secondary structure element) and those between sites in adjacent layers (forming the tertiary structure). The results, shown in Fig. 10, reveal that in all cases (the two contact distributions have again been combined in each plot) the contacts within layers tend to form first, followed, after a delay, by contacts between layers. Furthermore, the formation of contacts within layers is essentially complete, and the missing contacts are those that should have appeared between the layers. The detailed distributions, and especially the elapsed time between the formation of specific fractions of intra- and inter-layer contacts, vary substantially with configuration type (and, of course, with chain length).

The tendency for intra-layer contacts to form first, and in some instances for the core sites to form contacts earlier, are typical of the kinds of general observations that can be made concerning the folding pathways of simplified models. There is no component in the potential energy function that directly prefers any particular feature over any other (unlike the helix pairs of Rapaport (2002) where the potential favored the prior formation of individual helices and only subsequently their parallel alignment), so that this behavior is something that emerges spontaneously from the model. Owing to the wide variation in possible pathways, itself a consequence of the many degrees of freedom of the chain and the resulting high dimensionality of configuration space, it is not obvious how to further categorize the folding pathways at a higher level of resolution (by, for example, introducing appropriately defined “reaction coordinates”). In particular, it does not seem possible to establish the existence of “funnels” in configuration space that are traversed by a relatively large proportion of pathways; whether such a concept is useful (except in certain very specific cases) is questionable, indeed there could well be so many funnels that they span a substantial portion of configuration space.

### III.4 Contact maps

The contact map provides a highly detailed but concise summary of the proximity of a configuration to its native state by showing which site pairs are within contact range; it also reveals patterns in the distribution of contacts among chain sites as a function of their backbone separation and among groups of sites at particular backbone locations. It is also possible to forgo these details and study how the contact formation rate depends on backbone separation alone, without taking into account which sites are involved; such results reveal, not unexpectedly, that nearby (more “local”) sites generally form contacts earlier and have a higher overall success rate (there is the occasional exception). This classification, however, ignores features that have already been found to be significant, such as contacts within and between layers. Although the full contact map is not subject to this shortcoming, it does have the opposite problem of an excess of detail, and because sample sizes are smaller due to each pair of sites being treated separately, the results are subject to increased statistical noise.

The general arrangement of a contact map is that both axes correspond to the indices of the sites along the backbone; if a pair of sites $`(i,j)`$ are within contact range then a spot appears at the corresponding coordinates. Since the plot is symmetric only half the spots need be shown; the backbone pairs can be included for reference even though they are always present. It is possible to extend this scheme for representing the configuration to produce a reasonably reliable, low-temporal-resolution history of the folding process; this is accomplished by replacing the binary (pair/nonpair) values in the contact map by continuous values representing the fraction of configurations sampled in which pairing occurs (including only those sites that are actually paired in the native state) over the course of folding, and using either color or greyscale coding to denote the value. On the assumption that the contacts that appear more frequently are those that develop earlier along the folding pathway, the history is contained in this generalized version of the contact map; a considerably less concise alternative would be a sequence of binary maps corresponding to different times during the simulation.

The contact maps shown in Fig. 11 incorporate two different kinds of information. The first is the actual pattern of contacts appearing in the ground state of the chain, which is indicated by the presence of spots irrespective of their shading; only the results for the shorter chains are shown because for longer chains the details are too fine to be reproduced in these small figures. Certain structural features lead to recognizable patterns in the contact map, and while particular motifs appear in more than one map, the overall patterns associated with the four configurations are seen to be very different (there is also a systematic dependence on chain length). Contacts involving nearby (in terms of backbone separation) pairs appear close to the diagonal; identifiable structural elements in the native state, such as zigzag layers, spirals and reversed adjacent layers, correspond to specific spot patterns.

The other category of information present in the generalized contact map is a concise representation of the folding history of the chain ensemble. Since the greyscale representation embodies the averaged folding history, the time at which a contact typically appeared can be inferred from the darkness of the corresponding spot (the use of color increases the temporal resolution, so that color plots were also used in examining the behavior). What sort of information can be extracted from maps of this kind?

Given that in the present model the contacts themselves do not contribute directly to the interaction energy, but only reflect the organization dictated by the torsional potential, it is not surprising that the maps reveal that contacts between nearby sites along the backbone tend to appear earlier (there are, for example, no attractive forces acting directly between distant sites that could compete against this ordering). The behavior is not always monotonic, so that in the longer, linear spot sequences perpendicular to the main diagonal, a feature associated with configurations of type B and C, both of which have reversed adjacent layers, there are instances of more distant backbone sites pairing (on average) prior to those that are nearby. One could attempt to correlate specific features in the contact maps with the ability of the chains to fold successfully; while there are particular details that are likely to be relevant, such as the presence or absence of larger groups of spots at various distances from the main diagonal, the data are probably insufficient to allow any firm conclusions. What is apparent here is that the best folder, type A, has the least amount of structure of this kind, B and C are similar, while D, the poorest folder, has the most structure (these observations also apply to the longer chains). The inclusion of additional native-state structures in the study might help elucidate the relationship between contact patterns and foldability.

## IV Conclusions

The focus of this paper has been on the folding ability of model chain molecules with well-defined native states. While each of the native-state configurations is cube filling, the paths taken by the chains through the cubes are different for each configuration, and this difference – which in some sense determines the “accessibility” of the native state – has a strong influence on the outcome of the folding process. Unlike the requirement of many lattice studies, the native-state energy is not markedly lower than for compact misfolded states; nevertheless, some of the configurations already achieve high folding success rates under the conditions of the simulation, and the results for the less successful cases could be further improved by slower cooling. The key conclusion is, therefore, that if a unique native configuration exists that is consistent with the interaction potential, the chain is able – both in principle and in practice – to find it. The difficult problem, in the context of more realistic models, is how to formulate such a potential.

The relevance of the study of cubic conformations, despite the fact that such systems have no counterpart in nature, is the existence of well-defined, unambiguous ground states. Additionally, the cubic native state is able to accommodate a combination of both short- and long-range structural features. The fact that there have been extensive studies of cubically packed chains using other simulational approaches is yet another motivating factor. The advantage of the continuum MD approach is that efficient collective reorganization can occur, even in relatively compact states, as opposed to the discrete sets of moves provided by the various Monte Carlo approaches; the absence of configurational restrictions imposed by an underlying lattice allows for continuous conformational change, an important capability for partially-collapsed chains.

The most obvious limitations of the model, in its present form, is that the only interactions included, aside from excluded volume, are of torsional type, and that no solvent is included. Neither shortcoming is difficult to overcome in principle. Additional interactions could be included if there was a specific reason to do so, although consistency of the overall potential with the desired native state would have to be assured. The inclusion of an inert solvent, modeled using discrete particles, would slow the dynamics substantially and require additional computation to handle the extra degrees of freedom; it would also be possible to incorporate specific chain–solvent interactions (to mimic hydrophobicity for example), but this, too, would be a complicating factor that the deliberately simple design of the present model is aimed at avoiding.

Even with these simplification, the model has been shown to display a range of physically interesting and relevant kinds of behavior. All the observations are based on a chain model whose interactions have been designed without any specific folding scenario in mind (aside from the chosen native state); nevertheless, there are certain identifiable aspects of the mean behavior that merit attention, such as the tendency in certain configurations for the core sites to form contacts earlier, for layers to become organized first, and for the correct folding pathway to be correlated with a slower overall collapse (as measured by the radius of gyration). Such features correspond to generic behavior, so that simple models of this kind can be used as reference systems when investigating more highly refined designs. If additional features are able to produce little more than already occurs in the simplest of models, then they clearly fail to advance the state of the art.

What lessons emerge from these results insofar as protein folding is concerned? Each of the cases studied has a single minimum (free-) energy collapsed state and a variety of misfolded states with very similar energies. Thus there is no global minimum state that is well-separated from the many local minima, and on the basis of energy considerations alone it is not possible to reach any conclusions about the folding capability of each kind of chain. The intuitive but non-quantifiable notion of accessibility – the ability of the chain to go where it needs to be – plays an apparent role, since it reflects the robustness of the chain in resisting the consequences of incorrect entanglement; partially entangled states occur often along a folding pathway, but their effect is transient if escape is possible. It is clear from the present results that even the simple chain models studied here differ greatly in this respect.

###### Acknowledgements.

This work was partially supported by the Israel Science Foundation.


format: text

text:
# The Structure of Electroweak Corrections Due to Extended Gauge Symmetries

## 1 Introduction

Understanding the dynamics of mass generation requires looking beyond the standard model of particle physics. Models with extended electroweak gauge sectors play a large role in recent work on the origins of electroweak symmetry breaking, the large top quark mass and the top-bottom mass splitting . Some theories include an extended $`SU(2)\times SU(2)`$ structure for the weak interactions; the fermions may transform under the groups in a variety of flavor non-universal ways. Examples include non-commuting extended technicolor (NCETC) , topflavor , and the ununified standard model (UUM)  and related theories . Other theories include an extended $`U(1)\times U(1)`$ structure. Examples include topcolor-assisted technicolor , flavor-universal topcolor-assisted-technicolor , topflavor-seesaw , and string theories  or string-inspired models .

This paper studies the phenomenology of models with replicated electroweak gauge symmetries, based on a framework with the gauge structure $`SU(2)\times SU(2)\times U(1)\times [SU(2)\mathrm{or}U(1)]`$. We discuss flavor-universal theories of this type in detail and show how to extend our results in a straightforward manner to flavor non-universal theories.

We begin, in section 2, by establishing the conventions we shall use and demonstrate that many specific models discussed in the literature are special cases of our general framework. In section 3, we examine the corrections the extended gauge sector causes in the low-energy current-current interactions of both the charged and neutral gauge sectors. We find that the form of the corrections depends neatly on the charge assignments of the fermions under the extended electroweak group. In section 4, we perform a similar analysis of the masses and couplings of the $`W`$ and $`Z`$ bosons in our general framework. Again, the fermion transformation properties dictate the form of the leading corrections to the $`Zff`$ and $`Wff`$ couplings. In section 5, we perform a global fit of precision electroweak data to the theories in which the leading predicted corrections to observables are smallest. We discuss our conclusions and their application to other models in the literature in section 6.

## 2 Replicated Electroweak Gauge Symmetries

The general framework  for the models we will discuss is based on the gauge structure $`SU(2)\times SU(2)\times U(1)\times [SU(2)\mathrm{or}U(1)]`$, and is illustrated in the figure below. This figure employs the conventions of “moose” notation ,

(1)
in which the circles represent gauge groups with the specified gauge coupling, and the solid lines represent separate $`(SU(2)\times SU(2)/SU(2))`$ nonlinear sigma model fields which break the gauged or global symmetries to which they are attached. Each solid circle represents an $`SU(2)`$ or $`U(1)`$ group for which the entire global symmetry has been gauged; the dotted circle represents a global SU(2) in which only a U(1) subgroup has been gauged.

For convenience, the coupling constants of the gauge theories will be specified by

$$\stackrel{~}{g}^{}=\frac{e}{\mathrm{cos}\theta \mathrm{sin}\varphi },g^{}=\frac{e}{\mathrm{cos}\theta \mathrm{cos}\varphi },g=\frac{e}{\mathrm{sin}\theta \mathrm{cos}\omega },\stackrel{~}{g}=\frac{e}{\mathrm{sin}\theta \mathrm{sin}\omega },$$
(2)
and the $`f`$-constants (the analogs of $`f_\pi `$ in QCD) of the nonlinear sigma models by

$$\frac{f}{\mathrm{sin}\alpha },v,\frac{f}{\mathrm{cos}\alpha }.$$
(3)
The nonlinear sigma fields will break the extended gauge group first (at scale $`f`$) to an intermediate electroweak group $`SU(2)\times U(1)`$ and then (at scale $`v`$) to the unbroken electromagnetic group $`U(1)_{EM}`$. Hence the low-energy spectrum of the theory includes a massless photon, a set of $`W`$ and $`Z`$ bosons which closely resemble those of the standard model, and an extra set of heavy $`W`$ and $`Z`$ bosons corresponding to the generators broken at scale $`f`$.

As we will see, the Lagrangian parameters $`e`$, $`\theta `$, and $`v`$, will be approximately equal to the electric charge, weak mixing angle, and Higgs expectation value in the one-doublet standard model. The scale $`f`$ sets the masses of the extra gauge bosons, and the theory reduces to the standard model in the limit $`f\mathrm{}`$. The angle $`\alpha `$ allows us to independently vary the breaking scale of the duplicated $`SU(2)`$ or $`U(1)`$ gauge symmetries. The angles $`\varphi `$ and $`\omega `$ determine the couplings of the gauge bosons which become massive at scale $`f`$. This analysis will assume that the ratio $`\frac{v}{f}1`$ since precision electroweak data already tells us that the observed $`W`$ and $`Z`$ bosons closely resemble those of the standard model.

The symmetry structure of this model is similar to that proposed in the BESS (Breaking Electroweak Symmetry Strongly) model , an effective Lagrangian description motivated by strong electroweak symmetry breaking. The BESS model is in turn an application of “hidden local symmetry” to electroweak physics . Accordingly, we refer to our paradigm as “Generalized BESS” (GBESS). The symmetry structure in the limit $`f\mathrm{}`$ is precisely that expected in a “technicolor” model , and the theory has a custodial symmetry in the limit where $`g^{}`$ and $`\stackrel{~}{g}^{}`$ go to zero.

GBESS is the simplest model of an extended electroweak gauge symmetry incorporating both replicated $`SU(2)`$ and $`U(1)`$ gauge groups. As such the electroweak sector of a number of models in the literature can be thought of as forming special cases of GBESS. For example Noncommuting Extended Technicolor (NCETC) , can be framed in the moose notation as follows:

(4)
where the gauge groups correspond, from left to right, to the $`SU(2)_2\times SU(2)_1\times U(1)`$ electroweak structure of the model. In the “heavy” case of NCETC, the fermions transforming under $`SU(2)_2`$ are those of the first and second generations as shown, while the third generation transform under $`SU(2)_1`$; in the “light” case, the fermion assignments are reversed (per the double-headed arrow). The electroweak gauge structure of the Ununified Standard Model  is the same,

(5)
but it is the quarks and leptons which now transform under separate $`SU(2)`$ groups, reflecting the distinct weak interactions of these species in this model. As another example, the electroweak structure of topcolor , has the form

(6)
which illustrates the separate hypercharge interactions of the 3rd generation and light fermions.

Finally, electroweak $`SU(3)`$  has the gauge structure $`SU(3)\times SU(2)\times U(1)`$ which can be represented as follows

(7)
where the dashed circle reflects the fact that the gauged $`SU(2)\times U(1)`$ is a subgroup of a global $`SU(3)`$ group. Note that all fermions and the conventional Higgs field ($`\varphi `$) are charged under this $`SU(2)\times U(1)`$. The $`\mathrm{\Sigma }`$ transforms as a $`(2,3)_{\frac{1}{2}}`$ of the gauge symmetries. If we focus on the $`SU(2)\times U(1)`$ subgroup of the $`SU(3)`$ gauge group which is preserved by the $`SU(3)\times SU(3)/SU(3)`$ chiral symmetry breaking structure, this model can be redrawn.

(8)
This unfolds to reveal the linear moose form of GBESS once more, with the coupling values shown below

(9)
and with the standard model fermions transforming under the inner two groups. Note that for $`g,g^{}g_3`$ one finds $`\mathrm{sin}^2(\theta _W)\frac{1}{4}`$. The general properties of precision electroweak constraints on these models can correspondingly be viewed as special cases of what follows, as we will discuss below.

## 3 Low-Energy Interactions

Important constraints on models with extended electroweak symmetries arise from low-energy measurements, of which the most sensitive are the muon lifetime (used to determine $`G_F`$), atomic parity violation (APV), and neutrino-nucleon scattering. In the usual fashion, we may summarize the low-energy interactions in terms of four-fermion operators. The form of these interactions will depend, however, on the fermion charge assignments. For simplicity, we focus on models in which the fermion charge assignments are flavor universal. To illustrate the model-dependence of the results, we consider two examples.

First, we consider the case (hereafter, Case I) in which the ordinary fermions are charged only under the two groups at the middle of the moose

$$\text{}.$$
(10)
In this case the charged current interactions may be computed to be

$$L^{CC}=\frac{2}{v^2}J^{\mu +}J_\mu ^{},$$
(11)
and the neutral current interactions to be

$`L^{NC}=`$ $``$ $`{\displaystyle \frac{2}{v^2}}(J_3^\mu Q^\mu \mathrm{sin}^2\theta )^2`$ (12)
$``$ $`{\displaystyle \frac{2}{f^2}}(\mathrm{cos}^2\alpha \mathrm{sin}^4\theta \mathrm{sin}^4\omega +\mathrm{sin}^2\alpha \mathrm{cos}^4\theta \mathrm{sin}^4\varphi )Q^\mu Q_\mu .`$
In these expressions, the currents $`J_{\pm ,3}^\mu `$ and $`Q^\mu `$ are the conventional weak and electromagnetic currents. From these, we see that the strength of $`G_F`$, APV, and neutrino scattering is determined by $`v`$ in the usual way. Furthermore, comparing the two equations, we see that the strength of the charged and neutral current interactions, the so-called low-energy $`\rho `$ parameter, is precisely one (at tree-level). This last fact is a direct consequence of the Georgi-Weinberg neutral current theorem .

As an alternative, consider the case (hereafter, Case II) in which the $`SU(2)`$ charges of the ordinary fermions arise from transforming under the gauge group at the end of the moose ($`SU(2)_L`$ with coupling $`\stackrel{~}{g}`$)

$$\text{}.$$
(13)
A calculation of the charged current interactions yields

$$L^{CC}=2\left(\frac{1}{v^2}+\frac{\mathrm{cos}^2\alpha }{f^2}\right)J^{\mu +}J_\mu ^{},$$
(14)
while the neutral current interactions are summarized by

$`L^{NC}=`$ $``$ $`{\displaystyle \frac{2}{v^2}}(J_3^\mu Q^\mu \mathrm{sin}^2\theta )^2{\displaystyle \frac{2\mathrm{cos}^2\alpha }{f^2}}(J_3^\mu \mathrm{sin}^2\theta \mathrm{sin}^2\omega Q^\mu )^2`$ (15)
$``$ $`{\displaystyle \frac{2\mathrm{sin}^2\alpha }{f^2}}(\mathrm{cos}^4\theta \mathrm{sin}^4\varphi )Q^\mu Q_\mu .`$
Several points in this expression are of particular note. First, the value of $`G_F`$ as inferred from muon decay is no longer simply proportional to $`v^2`$, but is instead proportional to $`(v^2+\mathrm{cos}^2\alpha f^2)`$. As we shall see in the next section, this ultimately will give rise to corrections to electroweak observables of order $`(v/f)^2`$ which are not suppressed by any ratios of coupling constants. Second, in contrast to Case I, the strength of low-energy charged- and neutral-current interactions are no longer the same. The strengths of the $`J_3J_3`$ and $`J^+J^{}`$ portions of the interactions are still equal; this reflects the approximate custodial symmetry of the underlying model. The ratio of the strengths of the $`J_3Q`$ and $`J^+J^{}`$ terms, however, is no longer simply proportional to $`\mathrm{sin}^2\theta `$, but depends also on the angle $`\omega `$ relating the gauge couplings $`g`$ and $`\stackrel{~}{g}`$.

One can likewise define models in which the fermions couple to the group at the right-hand end of the moose instead of to the gauged $`U(1)`$ group on the interior. Whether the global $`SU(2)_R`$ is fully gauged or only a $`U(1)`$ subgroup is gauged, the effect on the neutral-current weak interactions will be the same; the charged weak bosons of a gauged $`SU(2)_R`$ would not mix with those of $`SU(2)_W`$ and $`SU(2)_L`$. We will refer to the situation where fermions are charged under groups W and R as Case III and that where fermions are charged under groups L and R as Case IV. The low-energy charged currents of Case III (Case IV) are identical to those of Case I (Case II). The low-energy neutral currents may be written in the general form:

$$L^{NC}=\frac{2}{v^2}(J_3^\mu Q^\mu \mathrm{sin}^2\theta )^2\frac{2\mathrm{cos}^2\alpha }{f^2}(t_3J_3^\mu t_QQ^\mu )^2\frac{2\mathrm{sin}^2\alpha }{f^2}(u_3J_3^\mu u_QQ^\mu )^2$$
(16)
where the coefficients $`t_3,t_Q,u_3,u_Q`$ are summarized in Table 1. Note that in Cases III and IV the strenghts of the $`J_3J_3`$ and $`J^+J^{}`$ portions of the interactions are no longer equal.

## 4 Z-Pole Constraints - General Structure

Many of the most significant constraints on physics beyond the standard model arise from precise measurements at the Z-pole. To interpret these measurements, we must compute the masses of the light $`W`$ and $`Z`$ bosons (which are the same in Cases I-IV) and their shifted couplings to ordinary fermions in terms of the Lagrangian parameters.

The charged weak bosons of the $`SU(2)_L`$ and $`SU(2)_W`$ groups of GBESS mix with one another222If one were to gauge the full global $`SU(2)_R`$ of GBESS, the charged weak bosons of group $`R`$ would not mix with those of groups $`L`$ and $`W`$ and would simply have masses $`M_{R^\pm }^2=\frac{e^2f^2}{4\mathrm{cos}^2\theta }\left(\frac{1}{\mathrm{sin}^2\alpha \mathrm{sin}^2\varphi }\right)`$. . Diagonalizing the mass-squared matrix

$$\frac{e^2}{4\mathrm{sin}^2\theta }\frac{f^2}{\mathrm{cos}^2\alpha \mathrm{sin}^2\omega }\left(\begin{array}{cc}1& \mathrm{tan}\omega \\ \mathrm{tan}\omega & \mathrm{tan}^2\omega \left[1+\frac{v^2\mathrm{cos}^2\alpha }{f^2}\right]\end{array}\right)$$
(17)
for the charged $`L`$ and $`W`$ bosons yields masses for the light $`W`$

$$M_W^2=\frac{e^2v^2}{4\mathrm{sin}^2\theta }\left(1\mathrm{cos}^2\alpha \mathrm{sin}^4\omega \frac{v^2}{f^2}\right)+O\left(\frac{v^4}{f^4}\right),$$
(18)
and its heavy counterpart (at leading order)

$$M_{L^\pm }^2=\frac{e^2f^2}{4\mathrm{sin}^2\theta }\left(\frac{1}{\mathrm{cos}^2\alpha \mathrm{sin}^2\omega \mathrm{cos}^2\omega }\right)$$
(19)
All four neutral bosons in GBESS mix with one another. Diagonalizing the mass-squared matrix (from which an overall factor of $`[e^2f^2/4\mathrm{cos}^2\theta \mathrm{sin}^2\theta ]`$ has been removed)

$$\left(\begin{array}{cccc}\frac{\mathrm{cos}^2\theta }{\mathrm{cos}^2\alpha \mathrm{sin}^2\omega }& \frac{\mathrm{cos}^2\theta }{\mathrm{cos}^2\alpha \mathrm{sin}\omega \mathrm{cos}\omega }& 0& 0\\ \frac{\mathrm{cos}^2\theta }{\mathrm{cos}^2\alpha \mathrm{sin}\omega \mathrm{cos}\omega }& \frac{\mathrm{cos}^2\theta }{\mathrm{cos}^2\omega }\left[\frac{v^2}{f^2}+\frac{1}{\mathrm{cos}^2\alpha }\right]& \frac{v^2}{f^2}\frac{\mathrm{cos}\theta \mathrm{sin}\theta }{\mathrm{cos}\omega \mathrm{cos}\varphi }& 0\\ 0& \frac{v^2}{f^2}\frac{\mathrm{cos}\theta \mathrm{sin}\theta }{\mathrm{cos}\omega \mathrm{cos}\varphi }& \frac{\mathrm{sin}^2\theta }{\mathrm{cos}^2\varphi }\left[\frac{v^2}{f^2}+\frac{1}{\mathrm{sin}^2\alpha }\right]& \frac{\mathrm{sin}^2\theta }{\mathrm{sin}^2\alpha \mathrm{sin}\varphi \mathrm{cos}\varphi }\\ 0& 0& \frac{\mathrm{sin}^2\theta }{\mathrm{sin}^2\alpha \mathrm{sin}\varphi \mathrm{cos}\varphi }& \frac{\mathrm{sin}^2\theta }{\mathrm{sin}^2\alpha \mathrm{sin}^2\varphi }\end{array}\right)$$
(20)
yields a massless photon, a light $`Z`$ state with mass

$$M_Z^2=\frac{e^2v^2}{4\mathrm{sin}^2\theta \mathrm{cos}^2\theta }\left(1(\mathrm{cos}^2\alpha \mathrm{sin}^4\omega +\mathrm{sin}^2\alpha \mathrm{sin}^4\varphi )\frac{v^2}{f^2}\right)+\mathrm{}$$
(21)
and two heavy $`Z`$ states with leading-order masses

$$M_{N1}^2=\frac{e^2f^2}{4\mathrm{cos}^2\theta }\left(\frac{1}{\mathrm{sin}^2\alpha \mathrm{sin}^2\varphi \mathrm{cos}^2\varphi }\right)$$
(22)
$$M_{N2}^2=\frac{e^2f^2}{4\mathrm{sin}^2\theta }\left(\frac{1}{\mathrm{cos}^2\alpha \mathrm{sin}^2\omega \mathrm{cos}^2\omega }\right)$$
(23)
By comparing our expression for $`M_Z^2`$ with our previous results for the relationships between $`G_F`$ and $`v`$ in Cases I and II, we can see that the fermion charge assignments have a major effect on the predicted values of electroweak observables. The leading corrections to the standard model relation between $`G_F`$, $`\alpha `$, $`M_Z^2`$, and the appropriately defined weak mixing angle $`\mathrm{sin}^2\theta _W`$ are generically of order $`v^2/f^2`$ in Case II. In Case I, the dominant corrections are further suppressed by factors of $`\mathrm{sin}^4\omega `$ or $`\mathrm{sin}^4\varphi `$. Then, if one views the electroweak predictions of GBESS in terms of corrections to the corresponding standard model results, these corrections in Case I are (potentially) suppressed by ratios of coupling constants relative to the size of corrections in Case II. This leads generically to weaker constraints in Case I models, meaning that new phenomena associated with such models (such as extra gauge bosons) can exist at lower energy scales.

In order to make predictions for electroweak observables, we need to compute the couplings of the ordinary fermions to the light gauge boson eigenstates. For Case I theories, we find that the couplings to the left-handed fermions to the $`W`$ are

$$\frac{e}{\mathrm{sin}\theta }\left(1\mathrm{cos}^2\alpha \mathrm{sin}^4\omega \frac{v^2}{f^2}\right)+\mathrm{}$$
(24)
and for the $`Z`$ we find the couplings

$`{\displaystyle \frac{e}{\mathrm{sin}\theta \mathrm{cos}\theta }}\left[1(\mathrm{cos}^2\alpha \mathrm{sin}^4\omega +\mathrm{sin}^2\alpha \mathrm{sin}^4\varphi ){\displaystyle \frac{v^2}{f^2}}+\mathrm{}\right]T_3`$
$`{\displaystyle \frac{e}{\mathrm{sin}\theta \mathrm{cos}\theta }}\left(\mathrm{sin}^2\theta \mathrm{sin}^2\alpha \mathrm{sin}^4\varphi {\displaystyle \frac{v^2}{f^2}}+\mathrm{}\right)Q.`$ (25)
Examining our expressions for the shifted gauge-boson masses and couplings, we see that, for Case I, all corrections to standard model predictions may be expressed in terms of two combinations of Lagrangian parameters:

$$c_1=\mathrm{cos}^2\alpha \mathrm{sin}^4\omega \frac{v^2}{f^2},c_2=\mathrm{sin}^2\alpha \mathrm{sin}^4\varphi \frac{v^2}{f^2}.$$
(26)
This allows us to compute bounds on model parameters in terms of fits to $`c_1`$ and $`c_2`$, greatly simplifying the calculations.

For Cases II-IV, the leading corrections to gauge-boson couplings at order $`v^2/f^2`$ are suppressed by only two powers of $`\mathrm{sin}\omega `$ or $`\mathrm{sin}\varphi `$. To facilitate comparison of the forms of the couplings in the different Cases, we write the general formula for $`W`$ couplings to left-handed fermions

$$\frac{e}{\mathrm{sin}\theta }\left(1+a_\pm \frac{v^2}{f^2}\right)$$
(27)
and for $`Z`$ couplings to fermions

$$\frac{e}{\mathrm{sin}\theta \mathrm{cos}\theta }\left[1+a_3\frac{v^2}{f^2}\right]T_3\frac{e}{\mathrm{sin}\theta \mathrm{cos}\theta }\left[\mathrm{sin}^2\theta +a_Q\frac{v^2}{f^2}\right]Q.$$
(28)
and tabulate the factors $`a_\pm ,a_3,`$ and $`a_Q`$ in Table 2. The general pattern is that when the fermions couple to a group on the interior (exterior) of the moose, the leading correction at order $`v^2/f^2`$ goes as $`\mathrm{sin}^4`$ ($`\mathrm{sin}^2`$) of the associated mixing angle. Thus the corrections in models where fermions couple to an exterior group can be significantly larger than those in Case I.

## 5 Flavor-Universal Analysis and Results

From the calculations in the previous two sections of the paper, we may compute the values of all precisely measured electroweak quantities in terms of the Lagrangian parameters. Using the procedure outlined in Burgess et. al., we find that the leading contributions to the observables are given in Table 3. As mentioned earlier, all correction terms of order $`\frac{v^2}{f^2}`$ are of the forms $`c_1`$ and $`c_2`$ defined in eqn. 26, meaning that they involve the fourth power of either $`\mathrm{sin}\omega `$ or $`\mathrm{sin}\varphi `$. In contrast, Table 2 illustrates the corrections of order $`\frac{v^2}{f^2}`$ for Case II-IV models are suppressed by fewer powers of these mixing angles; one expects terms proportional to $`\mathrm{sin}^2\omega `$ in Cases II, IV and to $`\mathrm{sin}^2\varphi `$ in Cases III, IV. We therefore perform a detailed comparison to the data only for models of Case I, in which all fermions are charged under the middle two groups of the moose.

We have made a global fit of the predictions of GBESS (Case I) to the electroweak observables listed in the most recent compilation by the LEP Electroweak Working Group . These include Z-pole observables as well as the width of the $`W`$ boson, and low-energy atomic parity violation and neutrino-nucleon scattering. The 68%, 95%, and 99% confidence level fits for the parameters $`c_{1,2}`$ are shown in Figure 1.

The best fit values are

$$c_1^{bestfit}=0.000563,c_2^{bestfit}=0.000108;$$
(29)
note that the negative value of $`c_2^{bestfit}`$ is unphysical given the definition 26 of $`c_2`$. The best-fit value of $`\chi ^2/df`$ is 30.5/20 for GBESS; the probability of a $`\chi ^2/df`$ of this size or greater is only 6.4%. Fitting the same data to the standard model, the $`\chi ^2/df`$ is 30.97/22, corresponding to a 9.7% probability. At present, no known model of this type gives a better fit to the data than the standard model.

Note that the constraints from the fit on $`c_2`$ (which depends on $`\mathrm{sin}\varphi `$) are significantly stronger than those on $`c_1`$ (which depends on $`\mathrm{sin}\omega `$). In a flavor-universal model, $`\mathrm{sin}^2\theta _{eff}`$ (the ratio of a fermion’s coupling to $`Q`$ to its coupling to $`T_3`$) is the same for each species. The majority of precision electroweak observables can be expressed in terms of ratios of the $`Z`$ boson’s couplings to various fermions – meaning that these observables depend only on $`\mathrm{sin}^2\theta _{eff}`$. Now the shift in $`\mathrm{sin}^2\theta _{eff}`$ from the standard model prediction (as inferred from $`G_F`$, $`\alpha `$, and $`M_Z`$ at tree level) depends as follows on $`\mathrm{sin}\varphi `$ and $`\mathrm{sin}\omega `$

$`\delta (\mathrm{sin}^2\theta _{eff})`$ $``$ $`\left({\displaystyle \frac{\mathrm{sin}^4\theta }{\mathrm{cos}^2\theta \mathrm{sin}^2\theta }}\right)\mathrm{cos}^2\alpha \mathrm{sin}^4\omega {\displaystyle \frac{v^2}{f^2}}`$ (30)
$``$ $`\left(1+{\displaystyle \frac{\mathrm{sin}^4\theta }{\mathrm{cos}^2\theta \mathrm{sin}^2\theta }}\right)\mathrm{sin}^2\alpha \mathrm{sin}^4\varphi {\displaystyle \frac{v^2}{f^2}}.`$
Since the first term in parentheses is $`0.1`$ and the second such term is $`1.1`$, then $`\mathrm{sin}^2\theta _{eff}`$ and the electroweak observables are far more sensitive to $`\mathrm{sin}\varphi `$ than to $`\mathrm{sin}\omega `$.

For a given value of $`\alpha `$, we may unfold these constraints to produce a 95% lower bound on $`f`$ in terms of $`\mathrm{sin}\omega `$ and $`\mathrm{sin}\varphi `$. A sense of the reach of these bounds is given in Figure 2, plotted for $`\alpha =\pi /4`$. For typical values of $`\mathrm{sin}\varphi `$ and $`\mathrm{sin}\omega `$, the bounds on the scale $`f`$ range from a few TeV. Note that the standard model structure is formally recovered in the limit where both $`\mathrm{sin}\varphi `$ and $`\mathrm{sin}\omega `$ approach zero.

## 6 Discussion and Conclusions

In studying the Generalized BESS class of models with extended electroweak gauge groups, we have seen how the assignment of fermion charges affects the structure of corrections to precision electroweak observables. Having fermions transform under the outer $`SU(2)`$ or $`U(1)`$ group of the linear moose (as in Cases II-IV) results in leading corrections that are suppressed by only two powers of the sine of the weak ($`\omega `$) or hypercharge ($`\varphi `$) mixing angle, respectively. If all fermions are charged only under the innermost groups (as in Case I), all corrections to electroweak observables are instead suppressed by $`\mathrm{sin}^4\omega `$ or $`\mathrm{sin}^4\varphi `$.

This fourth-power suppression was previously noted in the so-called degenerate BESS model  in which a parity symmetry imposed the condition $`\alpha =\pi /4`$. Our results show that the suppression is independent of $`\alpha `$, and arises purely from the choice of gauge charges for the fermions. Degenerate BESS exhibits the suppression not because of the parity symmetry (as originally supposed) but because its fermion couplings are as in Case I.

As noted earlier, a number of models in the literature are special cases of the GBESS framework. We can now understand the structure of the electroweak corrections found in the literature on these models in terms of the general properties outlined above. Many of these models are flavor non-universal, in contrast to the universal case discussed here. In such cases, one applies the rules outlined here separately to each species of fermion. For example, consider the heavy case of NCETC (equation 4). This is an example of GBESS in the limit where $`(\mathrm{sin}\alpha ,\mathrm{sin}\varphi )0`$ because there is no fourth group in the moose. The first and second generation fermions behave as in Case II of GBESS because they couple to the outer $`SU(2)_2`$ group; the 3rd generation fermions behave as in Case I. Therefore we predict the Fermi constant to be $`\sqrt{2}G_F=v^2+f^2`$ and the $`Z`$ coupling to left-handed fermions to be

$$\frac{e}{\mathrm{sin}\theta \mathrm{cos}\theta }\left[1+\frac{v^2}{f^2}(\mathrm{cos}^2\omega \mathrm{sin}^2\omega T_3^L\mathrm{sin}^4\omega T_3^H)\right]$$
(31)
where the separate weak generators for the heavy (H) and light (F) fermions are distinguished. Both of these results match what was found in . In the light case of NCETC, it is the light fermions which behave as in our Case I and the heavy fermions as in Case II; the predictions for $`G_F`$ and the $`Z`$ boson coupling to fermions are modified accordingly.333In making a comparison with light NCETC , one must note the following difference of conventions. In our discussion, we move fermions from one gauge group to another, leaving the order of the moose and the names of the couplings unaltered. In changing from heavy to light NCETC the original reference effectively moved the $`U(1)`$ hypercharge group to the left-hand end of the moose and swapped the couplings $`g`$ and $`\stackrel{~}{g}`$. Thus to compare with this reference we must make the exchange $`\mathrm{sin}\omega \mathrm{cos}\omega `$.

Many of the specific models cited above correspond to taking a weak-coupling limit for the extra gauge groups (i.e., taking $`\mathrm{sin}\varphi `$ or $`\mathrm{sin}\omega `$ of order 1), in which case the bounds on $`f`$ are of order 10 TeV, as consistent with results in the literature for topcolor, NCETC, UUM and electroweak $`SU(3)`$ . This tends to make the extra gauge bosons too heavy to observe easily in collider experiments. Formally the corrections to electroweak observables vanish when the couplings of the extra gauge groups become strong, that is in the limit $`(\mathrm{sin}\varphi ,\mathrm{sin}\omega )0`$. The phenomenologically interesting question is whether there are any viable models corresponding to this limit, which would result in interesting structure at accessible scales.

## Acknowledgments

Work partially supported by the US Department of Energy under grants DE-FG02-91ER40676 and DE-FG03-93ER40757. E.H.S. and R.S.C. acknowledge the hospitality of the Aspen Center for Physics where part of this work was completed.


text:
# Unexpected trends in the social facilitated survival of termites

## 1 Introduction

Grassé and Chauvin, (1944) reported a number of experiments to check the effects of group size on the survival of Apis sp., Polistes sp, Leptotorax sp., Formica sp. (Insecta: Hymenoptera) and Reticulitermes sp.(Insecta: Isoptera). In all of them, but in Polistes, survival increased with increasing group size. In Polistes they reported an “aberrant” result in which survival in both large and small groups was smaller than in the middle ones. The authors atributed these obeservations to bad experimental procedures or to a very particular physiological trait of the species. The experiment was hence repeated several times –even by independent observers– but the “aberrant” result was always found. On a survey experiment performed on the summer 1994, we test grouped individuals of Syntermes sp., Nasutitermes sp., and Cornitermes sp. (Insecta, Isoptera). Preliminary results showed the same unexpected behaviour on the survival of grouped Cornitermes cumulans (Kollar).

Such a convergence of results allow us to suspect that Grass & Chauvin (1944) have in fact unwarrantly spotted a pattern that is biologically sound. We base our suspicion on previous experimental works which show that (i) individual survivorship of termites, while affected by group size, may be related to individual mobility of group components (Miramontes and DeSouza,, 1996), and (ii) longevity of individual termites seems to obey a hump-shaped curve, clearly dependent on some (yet) unknown attribute of group size (DeSouza et al.,, 2001). In addition, there is theoretical evidence that cooperation could either stabilize or destabilize the dynamics of a social group, leading to maximal individual fitness at intermediate colony sizes (Avil s,, 1999).

This paper aims to present experimental evidence linking individual mobility, survivorship, and group size in termites. To achieve this we show that both individual mobility and survival are humped functions peaking at similiar group-size range. We argue, on experimental and theoretical grounds, that survival of starving grouped termite workers –at least in Cornitermes sp– obeys a humped-shape curve which present a biologically sound, not fortuitous, relationship to individual mobility.

## 2 Material & Methods

To confirm the trend observed in survey experiments, two further experiments were performed in 1996 and 1997 with grouped termites of this species collected from wild colonies in Vi osa in the state of Minas Gerais, Southeastern Brazil. The first experiment, spanning several days, was designed to measure longevity while the second was intended to measure individual activity. Both experiments involved termite workers (third instar and beyond) of different group sizes.

In the longevity experiment, termites from a wild colony were randomly placed in groups of 1, 2, 4, 8, 12, 16, 20, and 24 (each group replicated eight times). Collection took place in 04 April 1996. The groups were confined in test tubes made of transparent glass (9.5 $`\times `$ 1.4cm) with hermetically sealing rubber caps. Each test tube is a ‘replicate’. Test tubes have been previously washed and rinsed, soaked in sodium hypochlorin for 24h, rinsed again, and sterilised at $`180^{}\mathrm{C}`$ for one hour. Tubes with termites were kept horizontally, separated by plastic foam to prevent stridulation or other mechanically transmitted signals to propagate between the tubes. The workers were incubated in the dark at a constant temperature chamber ($`25^{}\mathrm{C}\pm 0.5`$) and were allowed to acclimatise for 12h. No food or water was provided, but tubes were opened each 24h, to allow air exchange. Tubes were exposed to light during the counting of survivors only (no more than 5 min). Observations were made each eight hours intervals, until all individuals were dead (153h; 20 observations). Termite groups in which dead individuals presented any sign of cannibalism were excluded from data analysis.

In the activity experiment, termites from another wild colony were randomly placed in groups of 2, 4, 6, 8, 16, 20, and 30 individuals (each group replicated four times). Collection took place in 12 October 1997. The experiment was setup as above, using a controlled temperature room ($`24^{}\mathrm{C}\pm 0.5`$) rather than a chamber. Termites were allowed to acclimatise in such a room, for three hours before the observations began. Two observers took notes of the number of termites moving in each group, at five minutes intervals, during ten continuous hours. Each observer managed two replicates of the experiment (two tubes of each group size). Light intensity in the room was the lowest possible to allow the work of the observers.

Statistical analyses aimed to check whether group size would affect (i) survival and (ii) mobility of termites. Either quantities were treated as the $`y`$-variable in separate regression analyses, always with group size as the $`x`$-variable. Survival was measured as the average number of hours spent until individual termites die in a given group, which was calculated for each replicate (a given test tube containing termites), using Weibull frequency distribution (Crawley,, 1993; Pinder III et al.,, 1978). Survival functions obey the general form:

$$\mathrm{ln}S(t)=\mu ^\alpha t^\alpha $$
(1)
where $`S(t)`$ is the proportion of individuals from the initial cohort (the initial group of termites) still alive at time $`t`$, $`\alpha `$ is the shape parameter, $`\mu `$ is the mean number of hours until death by the termites from a given replicate. The first step in the survival analysis was to establish the value of the shape parameter $`\alpha `$, thereby defining whether or not the proportion of individuals dying was constant through time ($`\alpha =1`$, $`dN/dt=1/\mu `$, exponential fit). To do so, survival curves were tentatively fit to the data, simulating different shape parameters, one at a time, until the error deviance reached an asymptote. The model obtained was then compared with a model assuming $`\alpha =1`$ (exponential fit), on the basis of their contribution to the error deviance. Parsimony requires to choose the exponential fit (simplest model) if it does not differ significantly ($`P`$ ¿ 0.05; $`\chi ^2`$ tables; 1 df) from an alternative, more complete, model.

After defining the shape parameter, the next step was to calculate the mean number of hours termites spent to die when confined in each group size ($`\mu `$ in equation 1 above). Each of such means represent the number of hours spent to die averaged across all individuals belonging to a given replicate, and are referred to as “mean time to death”. Mean times to death for each replicate were then collapsed into a single arithimetic mean for each group size, a valid procedure to avoid pseudoreplication effects (Crawley,, 1993). Such a procedure produces eight means (each one represented by a dot in Fig.1), one for each group size.

Means thereby produced were used to check whether termite survival (y-var) would be affected by group size (x-var). To do so, models were fitted to the data, starting from a null model ($`y=b`$, were $`b`$ is the grand mean), and adding new terms until achieving the best trade-off between percent of variance explained ($`r^2`$), and $`P`$-values obtained (see Table 1).

The general pattern of mobility of termites in a given replicate (a test tube containing termites) was measured by fitting a simple regression line of the form $`y=b+ax`$, through the data points formed by plotting the number of moving termites (y-var) against time in minutes (x-var). The slope $`a`$ of each line was used as a measurement of the general pattern of mobility for termites belonging to that replicate. If the number of moving termites increases as times goes by, the curve fitted to the data would show a positive slope. Conversely, a negative slope would indicate a decay in group mobility. When no clear trend is to be observed for a given replicate, its “mobility curve” would show a null slope. For each group size we obtained four slopes, since there were four replicates (i.e., four test tubes containing termites). The slope values thereby obtained were collapsed into a single arithimetic mean for each group size aiming, as above, to avoid pseudoreplication effects. Such a procedure produces seven means (each one represented by a dot in Fig.2), one for each group size. Such means are refereed to as the “general pattern of mobility” for that group size.

Means thereby produced were used to check whether the general pattern of mobility (slope $`a`$ of mobility lines, see above) would depend on the size of the group were termites have been confined. Models were fitted and selected as above (see Table 2).

All statistical analyses were performed using R (Ihaka and Gentleman,, 1996).

## 3 Results

The mean time to death of starved termites obeys a humped function of group size, presenting a peak at a characteristic group size and decreasing at both smaller and larger densities (Polynomial model, order four: F\[4;3\] = 20.22 ; $`r{}_{}{}^{2}adj`$ = 0.916; $`P`$ = 0.02; Fig. 1; Table 1) Accordingly, group activity presents a similar dependency on group size, showing a peak at an intermediate density and decreasing at both smaller and larger values (Quadratic model: F\[2;4\] = 9.095; $`r{}_{}{}^{2}adj`$ = 0.730; $`P`$ = 0.03; Fig. 2; Table 2). This quadratic model explains the data better than the simple linear one: adding a quadratic term to the linear model, increases by 18.7% the percentage of variance explained and improves the $`P`$-value (Table 2).

Such results are consistent across different procedures of data analysis. As specified above, each dot in Figs.1 and 2 represents an arithmetic mean taken across the respective replicates. If, however, we plot the raw data (i.e., splitting the arithmetic means into its components), we observe the same general shape in both curves (Figs.3 and 4).

## 4 Discussion

Across the Isoptera, a range of different socially facilitated behaviours has been extensively reported, although survival is the one that has attracted more attention (Table 3). In all survival studies, longevity has been reported to increase as group size increases. This may lead one to intuitively assume an asymptotic behaviour of such a curve, accepting that very large groups would attain the maximal possible survival rates. In this sense, the results reported by Grass & Chauvin(1944) in which survival of wasps is maximal at moderate densities (rather than at large group sizes), could be considered “aberrant”. However, a striking similarity can be observed between such results and those found here for C. cumulans termites (Fig. 1). More puzziling, termites from different wild colonies of the same species show a similar humped function, even when they have been poisoned by insecticide (DeSouza et al.,, 2001). Such a convergence of results lead us to suspect that this pattern is, in fact, biologically sound. It seems that survival of starving termite workers –at least in C.cumulans– is not asymptotically related to group size, but obeys a humped function, in which an optimal density (=number of individuals per unit area) assures maximal survival.

Optimal densities in social insects have already been predicted as a consequence of evolutionary pathways maximising individual fittness (Higashi and Yamamura,, 1993). Accordingly, extraordinary lifespan is argued to derive from evolutionary processes in which such insects benefit from sociality (Keller and Genoud,, 1977). Proximate mechanisms through which optimal densities would lead to better survival, however, are still poorly understood. Many works dealing with cluster size in group organisms are based on optimality arguments regarding resource usage, territory defence, anti-predatory behaviour and disease resistance (Creel,, 1997; Fritz and Garine-Wichatitsky,, 1996; Giraldeau and Beauchamp,, 1999; Roberts,, 1996; Rosengaus et al.,, 1998). Others, inspired on novel concepts of the sciences of complex systems, would invoke a relationship between group-size, information fluxes and task performance (Adler and Gordon,, 1992; Bonabeau et al.,, 1998; Deneubourg et al.,, 1986; Gordon,, 1996; Pacala et al.,, 1996). The interconnected nature of such systems is thought to be the basis for the self-organisation of a variety of cluster-related phenomen (Bonabeau et al.,, 1997; Camazine,, 1991; Miramontes et al.,, 1993; Theraulaz et al.,, 1995). Specifically for termites, evidences seem to point out that survival is related to mobility (Miramontes and DeSouza,, 1996). In addtion, studies addressing interindividual interactions and self-orgnisation predict the existence of group-size effects causing information transfer and behavioural diversity to reach near optimal conditions (Miramontes,, 1995; Solé and Miramontes,, 1995). Accordingly, termites studied here present maximal values for group activity at intermediate densities (Fig. 2), which lie in the same group-size range as the observed survival (Fig.1). It seems,therefore, that the association between group-size, non-asymptotic survival and greater individual activity is not fortuitous in C.cumulans. Moreover, to the extent that mobility affects the rate of social contacts –and thereby, cooperation– our results seem to agree to those of Avil s, (1999), who show an explicit link between cooperation, individual fitness, and group size in social organisms.

These results point out to specific mechanisms present in the dynamics of groups that may act to regulate observables such as density, number of individuals and rate of social contacts, that in turn have an impact on collective task performance. Arguably, behaviours such as colony fission, swarming, cannibalism, or resting, may be the natural mechanisms responsible for the regulation of the proportion of individuals in the colony. In other words, the patterns presented here may have unexpected consequences for understanding the life cycles of these organisms, and therefore are worth exploring further.

## 5 Acknowledgements

We thank CA Santos and DL Bernardo for assistance in laboratory and field work. Dr. A Chopps kindly refreshed our ideas with new and deep insights. OM would like to thank MP Hassel and P Rohani for assistance while visiting the UFViçosa and the Termitology Laboratory at UFV. ODS would like to thank the staff of the Departamento de Sistemas Complejos (IFisica, UNAM) for all the help during an academic visit. This work has been supported by UNAM-DGPA #INI08496, CONACyT #3280P-E9607, FAPEMIG, FINEP, and PROIN-Entomologia (CAPES-UFV). This paper was entirely produced using free software (Linux, , Xemacs, and R).

## 6 Legend to figures

Average number of hours termite workers spent to die when confined with conspecifics in test tubes, in the absence of food. Each dot represents an arithmetic mean across eight replicates (eight different test tubes). Average number of hours for each replicate was calculated by censored survival analysis, with Weibull distribution.

General pattern of mobility of termites confined with conspecifics in test tubes, in the absence of food. Mobility patterns are characterized here using the slope of a line fitted across a scatterplot of the number of moving termites (y-var) versus time (x-var), four each of the group sizes tested. Positive slopes denote that the number of moving termites increases with time, negative slopes denote the opposite. Null slopes denote a constant number of moving termites across the time range when observations were done. Each dot represents the average of the slopes presented by each of four replicates (four different test tubes).

Average number of hours termite workers spent to die when confined with conspecifics in test tubes, in the absence of food. Each dot represents the average number of hours termites survived in a single replicate (a test tube), as calculated by censored survival analysis, with Weibull distribution. Within each group size, datapoints of this figure were collapsed into a single arithmetic mean in order to produce Fig.1.

General pattern of mobility of termites confined with conspecifics in test tubes, in the absence of food. Mobility patterns are characterized here using the slope of a line fitted across a scatterplot of the number of moving termites (y-var) versus time (x-var), four each of the group sizes tested. Each dot represents the slope presented by a single of four replicates (four different test tubes). Within each group size, datapoints of this figure were collapsed into a single arithmetic mean in order to produce Fig.2.