text:
# Design of Toy Proteins Capable to Rearrange Conformations in a Mechanical Fashion

## I Introduction

### I.1 The problem and the context

Theoretical progress in understanding proteins in the recent years was concentrated on folding, along with connected questions of sequence design and evolution (see book fold\_book  and references therein for the recent overview). Folding attracts theorists not only because it is so important for fundamental biology and for pharmaceutical industry, but also because it is a robust universal phenomenon. Vast number of proteins exhibit the ability to fold, and there is a widely recognized necessity to understand the physical principle behind the selection of sequences capable of fast and reliable folding.

In our opinion, there is one more aspect of proteins which is equally robust and appealing for theoretical analysis in terms of some minimal model. We mean here the ability of many proteins to function in a machine-like fashion through ordered conformational rearrangement. This is most obvious for motor proteins whose function is directly related to certain mechanical (conformational) movements motor . This is also clear for ion channels whose function is to mechanically move molecules (or ions) from one place to the other (e.g., across the membrane). Although less obvious, conformational motions appear to be also very much at play in proteins whose function is purely electronic, such as, e.g., electron transfer in bioenergetics or catalysis of a chemical reaction. This latter point was first formulated by McClare McClare  and independently by Blumenfeld Blum\_Russian . More recently, it was extensively discussed in the book Blumenfeld . Somewhat different viewpoint on this subject was also recently presented in the book Chernavskii .

The new experimental data by H. Gruler et al Gruler  support the idea of slow conformational relaxation being important for the operation of enzimatic molecular machines. More detailed phenomenological models of enzyme operation based on the concept of conformational relaxation as a biased diffusion process have been successfully implemented to interpret the experimental results Mikhailov . There is now the Data Base of conformational movements in proteins database .

We emphasize two general properties of function-related conformational movements in proteins. First, they occur without significant opening of the dense globular structure. Viewed in the context of contemporary folding theories, this property seems quite exciting. Indeed, as native globule is pretty dense, frequently modeled as a maximally compact self-avoiding polymer DENSE , the inside movements may be expected to be strongly suppressed. Counterargument to this suggests that in fact real protein globule does have certain voids and is not absolutely dense VOIDS . Nevertheless, the density of a typical protein globule is similar to that of a polymer melt, for which reason it may be expected to be extremely viscous if not altogether glassy. The observation of significant conformational movements inside such a dense polymeric conglomerate challenges theory to offer an explanation.

Second general property we would like to mention is the presence of some preferred collective degree of freedom - which is almost a synonym to a functioning device. For instance, enzymes work in cycles, and each cycle means a turn around some loop in the conformational space. For channel-forming protein, a part of this loop corresponds to a transported ion moving from one place to the other, the rest of the loop corresponds to the protein coming back. According to the arguments developed a long time ago (see the book Blumenfeld  and references therein), it is important that there is only one collective degree of freedom along the loop of function (which, of course, does not rule out ”transverse” fluctuations fluctuations ).

Importantly, machine-like function is realized well away from equilibrium conditions. That means, there is no detailed balance, and the system moves along the loop in one direction and not in the other - there is no, and should not be, detailed balance. This, however, does not rule out the possibility that some parts of the cycle may present themselves as being the motions along the same path in opposite directions, like, e.g., a piston moving up and down in a steam engine. We shall return to this point later.

The notion of preferred function-related degree of freedom may be compared in some respects to the concept of reaction coordinate much discussed in folding studies trans\_coord ; DIRECT\_CURRENT . In both cases, the presence of transverse fluctuations is important. In case of folding, this gives rise to the understanding that, e.g., transition ”state” is not a microstate, determined to atomic details, but rather an ensemble of (micro)states Baldwin . The preferred functional degree of freedom must be considered in pretty much the same way.

It is an exciting question whether these two collective degrees of freedom, relevant for folding and function, are connected to one another. One could even speculate that they may be the same, or similar to some approximation. As of today, this question remains open. However, we note in passing that turn-around times reported in modern single molecule experiments on enzymes such as Feng\_Gai , on the order of a fraction of a millisecond, are not drastically different from typical folding times.

In the present paper, our goal is modest, but two-fold. First, we want to see, at least for the simplest model, how one can imagine a collective degree of freedom allowing orderly motion without opening the compact globule. Second, we want to design such a sequence that the globule energy changes in some desirable fashion while moving along the preferred coordinate. This way, we want to mimic a molecular spring.

Thus, our work is organized as follows. We first describe the model and formulate our problem in a more explicit way for that model. We then discuss the possibility to design the native state conformation, or, better to say, the ensemble of nearly-native state conformations, in such a way that they can realize the one-dimensional motion within an almost compact globule. After that, we design the sequence capable to fold into such a functioning conformation. At the end we study some properties of thus designed toy protein.

### I.2 A digression: from spherical cow to lattice protein

We shall work with the toy lattice model of protein. We understand that these words will run the emotions high and negative with many readers, and so we want to answer that from the very beginning. Everyone understands that there are no lattices in biological world, but this argument itself, however obvious, is too cheap to prove the lattice models useless. Indeed, for example, everyone laughs at the famous anecdote about a spherical cow, but at the same time everyone tacitly agrees that the model of a spherical cow is useful, e.g., to understand the scaling laws relating animal body mass and the rate of oxygen consumption NATURE . Thus, the dispute about usefulness, or the lack of one, for the lattice models in protein studies cannot be resolved on the level of philosophy, this is the question of specific purpose of certain studies. Not entering the details, there are questions for which lattice models are totally inappropriate (and may deserve laugh), there are some other questions for which using lattice models is legitimate. As we hope to prove by the results, our present paper belongs to the latter category.

Thus, we use standard toy lattice model in which protein is represented as a self-avoiding walk of the desirable length. The protein changes its conformation by means of elementary moves, including corner flips, end flips, crankshafts, and null moves VS . The advantage of this move set is that the resulting system is known to be ergodic ERGODICITY . For this model, polymer moves by making discrete succession of steps from one conformation to the next. Accordingly, the preferred degree of freedom must be associated with certain a linear (one-dimensional) succession of conformations, in which every conformation may only move into either previous or next conformation of the same group.

## II Design of conformation

### II.1 Space of conformations

We stated above why the protein needs to have a selected degree of freedom for functional work. Here, we design lattice toy protein in such a way that its conformations, while remain compact, can rearrange along one-dimensional linear path in conformational space.

The concept of linear path can be easily explained if the conformational space of lattice protein is visualized as a graph graph1 ; graph2 ; graph3 . In this representation, every conformation is denoted as a node of the graph, and two nodes are connected if and only if the transition between corresponding conformations is possible via single elementary move (see Figure 1).

We are looking for linear paths in the conformational space graph (CSG). Obviously, linear path is the succession of such nodes each of which is connected to exactly two other nodes. For swollen, non-compact polymer, a multitude of conformational motions is possible, the corresponding nodes of the graph have very many connections, and so none of the swollen conformations belong to the linear path. By constrast, in compact conformations conformational freedom of the polymer is very limited, and there is a hope to find linear paths.

In order to find one-dimensionally-connected paths of compact conformations, we examined from this point of view the properties of CSG of the short lattice polymer.

### II.2 Properties of the space of the compact conformations

We build on the findings of the work Foam . In that work, it was shown that placing the polymer inside a restricted size box makes the conformational space graph disconnected, consisting of several disjoint pieces, or chambers. This finding was based on computer simulation of lattice polymers of various lengths $`N`$, each confined in the $`3\times 3\times 3`$ box on the lattice. For our present purposes, it is important to address another physical situation, in which $`N`$ is fixed, while the degree of compactness of the chain may change. We achieve this by restricting the gyration radius of the chain $`R_g`$ and then looking at various specific values of $`R_g`$. More specifically, this was done in the following manner.

We consider lattice polymer of the length $`N=18`$. We start from maximally compact conformation and allow it to make all elementary moves consistent with the chain self-exclusion, but possibly (and necessarily!) violating the compactness. We accept the conformation and place it as a node on the graph if and only if the chain gyration radius in this conformation is less than the chosen threshold $`R_g`$. All the accepted conformations are pictured on the graph as nodes and their connections with all other accepted conformations are established through exhaustive search. Then these new conformations again allowed to make elementary moves, new conformations are accepted if they do not exceed the same threshold $`R_g`$, and the process repeats. As regards the limiting value of $`R_g`$, we choose it experimentally, and it regulates both the compactness of the conformations and the number of conformations in the graph. The graph constructed in such a manner consists of several disjoint regions Foam , or chambers. That means, for two conformations, which belong to different chambers, there is no sequence of elementary moves, which transforms one of them into another without breaking the restriction on $`R_g`$. Thus, the procedure must be repeated starting from different maximally compact conformations to list all the chambers of the graph. The number of conformations in different chambers varies, reflecting the distribution of the clusters in the bond percolation problem (in this case, we deal with percolation in conformation space Foam ). Figure 2 shows the dependence of the number of the small chambers in the conformational space graph on the number of conformations locked in chamber. In terms of the underlying physical idea, this figure is similar to the result reported earlier in the work Foam  for a different model. In the present model, we vary $`R_g`$ at the fixed number of monomers $`N`$. In the work Foam , the similar alleged percolation in the conformational space was controlled by changing $`N`$ while locking the polymer inside the $`3\times 3\times 3`$ cube on the lattice, which of course implies fixed $`R_g`$.

Further, we established the connectivities of the nodes which belong to the largest chamber in the CSG. The distribution of the connectivities of the nodes is shown in Figure 3a, curve 1. It is compared with the distribution of the numbers of neighboring nodes for the same set of conformations, but not restricted with the $`R_g`$-condition (curve 2). Curve 2 is indistinguishable from the binomial distribution (which does not contradict the idea that non-restricted CSG is a small-world network smallworld ). The CSG built under the restriction on the values of $`R_g`$ of toy-protein conformations (curve 1) is significantly different, one can imagine this graph as a percolation cluster of the bond percolation problem on the lattice with the topology of the small-world network percolation . For the small values of $`R_g`$ (weakly connected cluster) the peak of the distribution corresponds to the graph nodes connected to only two neighbors. The sharp peak on the curve 1 corresponding to the poorly connected conformations can be easily explained. The change of the geometry of the voids in the bulk of the protein globule is possible only via small number of local moves, because the excluded volume effect is very strong in compact conformations. Such subtle conformational moves do not affect significantly the value of the $`R_g`$ of the protein chain, whereas opening of some loop on the surface of the globule leads to the increase of the $`R_g`$. On the other hand, the majority of the conformational moves accessible to the unrestricted protein chain occurs on the surface. Accordingly, the limiting of $`R_g`$ from above forbids surface moves but does not restrict the changes in the bulk of the globule. That is why for small limiting $`R_g`$ values the sharp peak of the distribution rises at the poorly connected conformations.

Thus, the study of $`R_g`$-restricted conformations of the $`18`$-mer provides an example in which there is the peak of connectivity distribution which corresponds to the graph nodes with connectivity $`2`$. These are desired conformations forming linear paths. However, most of these linear paths are rather short and consist of only few elementary moves. The distribution of the lengths of the paths is shown in Figure 3b, it is exponential distribution.

Thus, our conclusion so far is this. Long linear paths exist, but they are not very common, they are rare. In this sense, the situation is similar to that of selection of sequences for toy proteins capable of folding fold\_book . In sequence selection case, the ”good” sequences are exponentially rare among the random ones. The goal of sequence design, or selection, is to fish them out. Similarly, our goal now is to identify the rare conformations which are connected by long linear paths in the conformation space. In order to do that, we need to understand better the local geometry of conformations belonging to the linear paths. It cannot be done for the chain length of $`N=18`$ which is too small. It is also small compared with typical protein lengths, about one hundred or more in average. Thus, we repeated the same procedure of CSG mapping for the toy-protein of the length $`N=125`$. Of course, in this case no exhaustive enumeration of conformations is possible, and so we performed random sampling instead.

We started from the limiting value $`R_g=2.44949`$, which is equal to the gyration radius of the maximally compact conformation. In this case, no conformational movement is possible, and, therefore, CSG consists of as many disconnected nodes as the number of compact conformations. We now increase $`R_g`$ by a very small amount. By several attempts, we choose the step of $`R_g`$ increment in such a way that after one step conformation space remains disintegrated, non-ergodic, consisting of disconnected chambers Foam  most of which contain only few conformations. We then increase $`R_g`$ step by step, every time mapping the newly obtained parts of conformation space on the graph, as described before. When $`R_g`$ approaches $`2.46`$, the percolation transition occurs and there appears an ”infinite” cluster. We mapped into the graph $`1000000`$ conformations from this cluster and searched for the long linear paths on the graph, as it was done in case of $`18`$-mer.

The longest paths we found for the $`125`$-mer are organized as follows. The conformations along the path transform into one another in a set of subsequent corner flips in the bulk of compact conformation. When the first corner flip occurs, the vacant site opens producing the opportunity for another monomer to move. Then the similar process repeats many times, leading effectively to the vacancy traveling through the globule, finding the new corner flips on its way. These subsequent transformations are demonstrated in Figure 4 which for the simplicity of presentation shows a smaller polymer, namely $`63`$-mer.

### II.3 Design of the pathway in conformation space

Now, when the mechanism of rearrangements within the compact conformation for the toy protein is clarified, we can propose the effective and straightforward approach to design the linear path of compact conformations.

The algorithm includes two main stages. First, we arrange the switching elements along the path. Second, the rest of the polymer is computationally designed such that it is (almost) maximally compact and contains the necessary set of switching elements. As a result, we obtain one-dimensional path of compact conformations of the lattice toy protein.

#### II.3.1 Building the switching elements.

We build the linear path using flipping corners as switching elements.

In the beginning of the procedure we have an empty cubic lattice. We start from choosing initial position of the vacancy (Figure 5a, node $`(0,0,0)`$). The first flipping corner is drawn on the lattice in such a way that when it flips, the corner and the vacancy exchange their positions. The edges forming next switching element should be drawn in such a way, that when it will flip, vacancy will hop to the former corner position opening the way for the next switching element. In the Figure 5b two subsequent flipping corners are shown. After they both will move, the vacancy will take the position $`(1,2,1)`$. Other flipping corners are drawn subsequently as many times as needed (Figure 5c).

At the same time as the switching elements arranged on the lattice we can control the linearity of the path. Every time when the new flipping corner is set on the lattice there are several ways to locate it relative to the vacancy position. We choose only one location, but the problem is that the other flipping corners may appear later, on the stage of the design of the whole conformation (see Figure 5d). To prevent appearance of such parasitic switches we draw additional pieces of the conformation surrounding the pathway sites, as it is shown in Figure 5e.

#### II.3.2 Building the rest of conformation.

We now have the set of switching elements, which are supposed to be just the disconnected pieces of the polymer. We have to find now the rest of the polymer, such that it fills (almost) completely the whole volume of the cube, and connects all switching elements into a linear chain. For this purpose, we developed computational method which is the modification of the approach proposed by Ramakrishnan et al Conf-des  to generate maximally compact conformations. Here, we describe our modified algorithm.

The conformation design starts with placing on the lattice the new edges connecting randomly chosen neighboring vertices (Figure 6 a,b). This process soon brings us to the state, where some vertex cannot be connected to its randomly chosen neighbor, because the neighbor already has two edges incident on it (Figure 6 c). In this case, special procedure called *two-matching* is applied Conf-des . During this and the following steps of the algorithm, some randomly chosen edges can be removed from the lattice and changed by others. However, we impose the condition that the edges forming switching elements *cannot be removed* at any step of conformation design.

*Two-matching* starts from picking up the vertex $`P`$, which is either not connected, or has only one incident edge. Then its neighbor $`Q`$ is chosen randomly as an opposite end of the new edge. If $`Q`$ belongs to the linear subchain, then the ends of this subchain are checked for the possibility to be connected with their neighbors. If it can be done, a new edge (the edge $`ST`$ at the figure 6b) is added. One of the edges incident on $`Q`$ is replaced by the edge $`PQ`$. Thus, in this procedure, the number of edges increases by one (Figure 6d). If, on the other hand, the vertex $`Q`$ belongs to the looped subchain, one of its incident edges is removed and replaced by $`PQ`$. So, the loop is broken, and the total number of edges on the lattice remains unchanged. Typically, as a result of the work of this procedure, we obtain several looped and one linear subchain packed into the cubic lattice (Figure 6e).

Now, the subchains should be merged into one chain. This is achieved in the following way. Suppose four neighboring vertices $`(K,L,M,N)`$ form a square. The connecting edges $`(K,L)`$ and $`(M,N)`$ belong to the different subchains. Excluding these edges and including instead $`(K,M)`$ and $`(L,N)`$, we would have merged subchains (see Figure 6f). Such an operation is known as  *patching*. Each patching operation transforms a pair of subchains into one subchain. The edges to be involved in patching are chosen randomly. The process is stopped when there is no more and no less than one linear chain on the lattice, which is the desired polymer conformation.

In the original work Conf-des , this method was applied to generate maximally compact lattice conformations. It is worth repeating that for our purposes we use this method starting from a complicated lattice which is the cube minus elements chosen for switching elements.

Generally, it is possible to use also crankshafts and flipping ends to design the switching elements. One just needs to forbid all the states of these elementary moves but two. This should be done by placing on the lattice additional edges (as it was done to prevent parasitic corner flips), which restrict the extra states of these moves. However, the end flip can be used only twice as switching element, because there are two ends of the chain. As regards the crankshaft move, it needs two vacant sites to make switching possible, which means the conformation in question should be slightly less compact. For these reasons, we use only corner flips as switching elements in this work.

## III Sequence design

After we have chosen conformations, we should find the sequence which fits the target conformations with low energies. For this purpose the sequence of the model protein is annealed and Monte Carlo optimization in the sequence space is performed. The details of the algorithm are previously published seq-des  (see also review article Grosberg1  and references therein for further details). It must be emphasized that we plan to work with the sufficiently large set of monomer species; in fact, we shall even use the so-called Independent Interaction Model, in which the number of distinct species is as large as the number of monomers in the chain. This allows us to avoid difficulties well known in the case of sequence design for the two-letter heteropolymers, such as the $`HP`$-model Yue . In the context of the present work, sequence design method had to be modified in two respects. First of all, we need not only one target conformation to have a low energy, as it is typically assumed in protein folding simulations. We need the whole family of conformations - all conformations belonging to the rearrangement path - to have distinctly lower energies than all other states. Second, the more ambitious goal is to design the sequence in such a way that moving of the system along its pathway changes energy in an orderly fashion. Since our rearrangement path is one-dimensional, we can pretend energy to increase monotonically along this path, in this sense making our toy protein a model of a molecular spring.

Let us consider these two aspects of sequence design one by one.

### III.1 Sequence with multiple ground states

The model protein is determined by the set of the coordinates of the monomers $`𝒞=\{\stackrel{}{r}_I\}`$ and the sequence of monomers $`seq=\{s_I\}`$ (the species $`s_I\{1,2,\mathrm{}q\}`$ denote the identity of each monomer, $`q`$ is the total number of species), index $`I`$ counts monomers along the chain. The Hamiltonian is written as follows

$$(\mathrm{seq},𝒞)=\underset{I<J}{\overset{N}{}}B_{s_Is_J}\mathrm{\Delta }(\stackrel{}{r}_I\stackrel{}{r}_J),$$
(1)
where the energy of the conformation is determined by the matrix of species-species energies $`B_{i,j}`$ for the contacting monomers and function $`\mathrm{\Delta }(\stackrel{}{r}_I\stackrel{}{r}_J)`$ is defined such that it is equal to $`1`$ if monomers $`I,J`$ are lattice neighbors and $`0`$ otherwise. We use independent interactions model (IIM) IIM  and Miyazawa-Jernigan (MJ) MJ  matrices for Monte Carlo simulations in this article.

In our approach the goal is to design sequences for which the whole set of conformations $`\{𝒞_k\}`$ have energies sufficiently below that of the rest of conformation space. Of course, our candidates for the target states $`\{𝒞_k\}`$ are the conformations which belong to the previously designed linear rearrangement path in conformation space. These conformations are supposed to form a deep valley in the energy landscape. For the set of two target conformations, sequence design was performed in the paper 2states . In that work, the goal was to model proteins which can fold into two (or more) distinct ”native” conformations, like prions Prusiner . Accordingly, two target conformations were chosen to be totally dissimilar (non-overlapping, or weakly overlapping). In more details, such design was examined in Multiple\_Funnels . In our case, the problem is almost the opposite. The conformations in question are very closely related, they can be mutually transformed into one another in just a few moves. Accordingly, the overlap between neighboring conformations along the path is very high. Of course, as the system walks along its pathway from one end to another, the overlap decreases, but still remains significant. For example, the overlap between the conformations at the opposite ends of the path shown in the Figure 4 is as high as about $`75\%`$.

The sequence optimization is governed by the following Hamiltonian:

$$_{\mathrm{des}}(\mathrm{seq},\{𝒞_k\})=\underset{k=1}{\overset{N_𝒞}{}}(\mathrm{seq},𝒞_k),$$
(2)
where $`N_𝒞`$ is the total number of conformations along the rearrangement pathway.

The question which arises now is this. How efficient is the sequence optimization in the case of multiple closely related target states? Let us consider the simple situation when only energies of two target conformations are optimized. These conformations, namely $`𝒞_1`$ and $`𝒞_2`$, could be, for example, the ends of the rearrangement pathway. How deeply can their energies be lowered during the sequence design in comparison with an arbitrary other compact conformation $`𝒞_{}`$?

To make this estimate we can calculate the energy of the conformation $`𝒞_{}`$ averaged over the sequence space

$$E_{}(seq)=\frac{_{seq}P_{seq}^{(0)}e^{\left(_{\mathrm{des}}(\mathrm{seq},\{𝒞_k\})\right)/T_{\mathrm{des}}}(\mathrm{seq},𝒞_{})}{_{seq}P_{seq}^{(0)}e^{\left(_{\mathrm{des}}(\mathrm{seq},\{𝒞_k\})\right)/T_{\mathrm{des}}}},$$
(3)
where $`P_{seq}^{(0)}=_{I=1}^Np_{s_I}`$ is the probability for the sequences made randomly from independent monomer species with occurrence probabilities $`p_i`$. The details of similar calculations are described in the review article Grosberg1 . The result for the present case of two target conformations reads

$$E_{}(seq)=\overline{B}Q^{}\frac{\delta \widehat{B}^2}{T_{\mathrm{des}}}[Q^{1,}+Q^{2,}],$$
(4)
where

$`\overline{B}={\displaystyle \underset{ij}{}}p_iB_{ij}p_j\mathrm{and}`$ (5)
$`\delta \widehat{B}^2={\displaystyle \underset{ij}{}}p_i(B_{ij}\overline{B})^2p_j`$ (6)
are the mean value and the variance of the interaction matrix,

$$Q^{}=\underset{I<J}{\overset{N}{}}\mathrm{\Delta }(\stackrel{}{r}_I^{}\stackrel{}{r}_J^{})$$
(7)
is the total number of contacts in the conformation $`𝒞_{}`$ and its overlap with arbitrary target conformation $`𝒞_k`$ is defined as

$$Q^{k,}=\underset{I<J}{\overset{N}{}}\mathrm{\Delta }(\stackrel{}{r}_I^k\stackrel{}{r}_J^k)\mathrm{\Delta }(\stackrel{}{r}_I^{}\stackrel{}{r}_J^{}).$$
(8)
As we can see from the expression (4), the energy of the designed sequence in the conformation $`𝒞_{}`$ depends on the similarity of this conformation to the target states. This similarity is measured by the overlap parameter $`[Q^{1,}+Q^{2,}]`$. It takes the maximal value for the a given pair $`𝒞_1`$, $`𝒞_2`$, when $`𝒞_{}`$ coincides with either $`𝒞_1`$ or $`𝒞_2`$. Usually, if $`𝒞_{}`$ lies on the rearrangement pathway, the overlap parameter $`[Q^{1,}+Q^{2,}]`$ takes values close to the maximum. Therefore, not only energies of the target conformations are optimized, but conformations between them (e.g. along the designed linear path!) are optimized, too.

In principle, sequence design method may also lower the energies of other states, which are related to the target states, but do not belong to the rearrangement pathway, and thus are non desirable for us here. This is a well known problem, generally addressed through the ”negative design” (see, for instance, negative ). Luckily, in our specific case, conformation design, as discussed above, helps to address this problem. Indeed, since there are no allowed conformations at all on the sides of the pathway, the only possible low energy decoys are structurally unrelated ones.

Sequence design method employed here may seem to contradict the results of the recent work manystates . In this work, authors estimated the maximal possible number of ”native” states which may be ”memorized” by the sequence. They showed that this number is very limited, it is independent on protein length, and is fully determined by the alphabet - the number, $`q`$, of distinct monomer species (for instance, for the system with $`q=20`$ monomer species, there can be not more than four or five ”native” states). In fact, there is no contradiction. The estimate of the work manystates  determines the maximal number of unrelated conformations which can be designed into the sequence. In our case, all target states are very closely related, and they in fact belong to the same potential well, or funnel, in the free energy landscape. They, of course, cannot be considered as independent. As we show below, our sequence design process works successfully even when the number of states in the rearrangement path is as large as about ten.

### III.2 Sequence design with energy gradient along the rearrangement path

So far, we have been discussing the sequence design procedure for which all the target states $`𝒞_k`$ were equal. Now, following Orwell Orwell , we shall consider some of the conformations more equal than others. Specifically, we should remember that conformations $`𝒞_k`$ form a one-dimensional path. We did it on purpose, and we should use it now. To be specific, let us assume that conformations $`𝒞_k`$ are labeled with index $`k`$ in a natural way, such that $`k`$ changes orderly from $`1`$ at the one end of the rearrangement pathway to the maximal value $`N_𝒞`$ at the other end. Then, we shall require that, say, $`𝒞_1`$ has the lowest energy, $`𝒞_2`$ has energy a little higher - preferably, by a certain amount $`W`$ higher; $`𝒞_3`$ we want to be about $`W`$ higher in energy than $`𝒞_2`$, … , and this continues all the way up to $`𝒞_{N_𝒞}`$, which we want to be about $`\left(N_𝒞1\right)W`$ above $`𝒞_1`$, but still much lower than all the other conformations.

Why do we want such energy landscape? First of all, being all connected in one valley, the conformations $`𝒞_k`$ form *together* a basin of attraction for folding, or folding funnel. Second, when correctly folded and at the bottom of the funnel, the system can still travel back and forth along the one-dimensional rearrangement path. This travel may be either due to fluctuations or due to some externally applied force.

To achieve this aim we modified the design Hamiltonian (2) in the following way:

$`_{\mathrm{des}}(\mathrm{seq},\{𝒞_k\})={\displaystyle \underset{k=1}{\overset{N_𝒞}{}}}(\mathrm{seq},𝒞_k)+`$ (9)
$`+`$ $`\lambda {\displaystyle \underset{k=2}{\overset{N_𝒞}{}}}\left[(\mathrm{seq},𝒞_k)(\mathrm{seq},𝒞_1)kW\right]^2,`$
where $`\lambda `$ is the ”experimentally” adjusted parameter, and $`W`$ is the desired energy gap between neighboring conformations. The Monte-Carlo optimization ruled by the Hamiltonian (9) has a bias towards sequences with lower energy in conformation $`𝒞_1`$ and subsequently higher energies in conformations $`𝒞_2,𝒞_3,\mathrm{}`$.

Various regimes are possible here depending on the relation between the design temperature $`T_{\mathrm{des}}`$, real temperature $`T`$, and the value of $`W`$. Not entering the discussion of all these regimes, we mention that in what follows we have chosen the sequence optimization temperature to be $`T_{\mathrm{des}}=0.1\delta B`$.

## IV Results

### IV.1 Design of the conformations and the sequence

To demonstrate the work of our design method, we generated two sets of conformations of lattice $`63`$-mer as described above in the section II.3. The linear rearrangement path of the conformation shown in Figure 4 includes 7 conformations. In another example, shown below in Figure 10, there are 6 conformations in the path. The location of the switching elements in these two cases is chosen differently. Switching elements in the former case are located in the bulk of the globule, whereas in the latter case the switching elements are all located near the surface.

For these two sets of target conformations, the sequence optimization procedure was applied. The values of parameters were as follows: $`W=0.3\delta B`$, $`T_{\mathrm{des}}=0.1\delta B`$, $`\lambda =15`$, and the interaction matrix was chosen to be correspond to the independent interactions model (IIM).

### IV.2 Folding

First of all, we have to check that the chains can correctly fold into the conformation as designed. We used the set of conformations shown in the Figure 44 as target states. We compared folding rates for the sequence designed as proposed in this article and for the control sequence, designed in a more traditional way, with conformation $`𝒞_1`$ as the only purported ground state. All the folding experiments were started from different random unfolded conformation. The Monte Carlo simulations were performed at temperatures in the range $`T=0.7÷0.9T_m`$, where $`T_m1.3\delta B`$ was the midpoint temperature. Mean first passage time (MFPT) at every temperature was calculated by averaging over $`30`$ folding runs. The results are shown in the Figure 7. The folding times for the sequence which has multiple ground states are approximately $`3`$ times longer than for the sequence with the unique native conformation. Further inspection suggests that this happens because the depth of the global minimum for the sequence with multiple ground states can not be as well optimized as that for the sequence with unique ground state. Nevertheless, the emphasis of our result here is on the good news, not the bad ones: it is not important that our sequences are slower, it is important that they are insignificantly slower, *only* by a factor of about $`3`$ slower. That means, they do fold, and their folding time is of the same order of magnitude.

It is worth to emphasize that during this folding experiments, the chain was not confined in any restricted volume. We used volume restriction in the preliminary stage of this work, to elucidate the method of conformation design. Now, as we are done with the design, we let the polymer do whatever it wants, and the result is that it folds and spontaneously arrives into the valley where it has a linear chain of conformations at its disposal.

### IV.3 Diffusion along the designed path

Now, since we have established that our model does fold, we have to check if it can move along the designed path. In reality, conformational relaxation of a functioning protein machine is triggered by the attachment or detachment of a substrate or other ligand; we shall consider this in the next section. Here, we want to perform the simpler test to see what happens without any stimuli. In this situation, we expect our toy protein to move randomly back and forth along the designed path. It should be mostly in the lowest energy state $`𝒞_1`$, but since the energy difference between states along the path $`W=0.3\delta B`$ is only a fraction of thermal energy, bias towards the $`𝒞_1`$ end should be relatively weak, there should be plenty of fluctuations along the path. The point to be checked is that the system, while performing this random walk along the path, should not open up too frequently, the globule should stay compact.

To examine how this happens in the toy-protein shown in the Figure 4, we run a long Monte Carlo simulations at different temperatures starting from the conformation $`𝒞_1`$. The events of passing the conformations $`\{𝒞_k\}`$ of the pre-designed path are recorded. Coordinate along the pathway takes the value $`x=k`$ if the vacancy position along the path coincides with conformation $`𝒞_k`$. A typical ”trajectory” of conformational changes for the first designed toy protein is shown in Figure 8. The inset of Figure 8 displays the details of one particular passage along the path from conformation $`𝒞_1`$ to $`𝒞_7`$ and back. The events when conformation changes in a way other than walking along the path (say, some loop opens on the surface of the globule), are pictured as the change of the coordinate in the perpendicular plane. As one can see, below the midpoint temperature the toy-protein stays steadily on the rearrangement path and makes random walks back and forth along it, with very limited fluctuations in the multitude of transverse directions in the conformation space.

Thus, what we observe confirms that our toy protein performs thermal fluctuations in the form of one-dimensional movement along the designed path. This is important, because if fluctuations occur that way, that means, under some different circumstances, e.g., applied external force, the system can perform also a forced movement along that same pathway, as it is expected for the function. Thus, we need to test carefully that the observed fluctuations are indeed the random walks effectively in one dimension. To address this quantitatively we calculate the correlation function $`x(t)x(t+\tau )`$ for the toy protein and compare it with that for the artificial auxiliary random walker which diffuses on the $`1D`$ discrete set of states with the same energy profile as for the protein. The good agreement between the two correlation functions is evident in Figure 9. Hence the designed toy protein does indeed move along the one-dimensional pathway and in this sense it remembers its ”function.”

### IV.4 Simulation of the binding/dissociation with the small molecule

In the previous section we demonstrated that toy-protein can use the designed rearrangement path, performing a biased, but random motion. This biased walk can be triggered if the protein molecule is externally stimulated. As in the real biochemical world, this can be most efficiently done by attaching the ligand molecule to the protein globule. In this sense, the typical cycle of a protein work may be described as follows Blumenfeld .

We start from the protein molecule in its ground state. At some moment, a ligand molecule gets attached to the protein. With the ligand attached, protein conformation is no longer the ground state, some other conformation is now lower in energy, and, therefore, attachment of ligand initiates the relaxation process of the globule into the new conformational state, which is the ground state for the protein-plus-ligand complex. When this relaxation is completed (or nearly completed), the new conformation turns out well suitable for certain chemical (or other small length scale) changes, the result of which for the protein is the desorbtion of a ligand initiating again the relaxation process, this time back into the original ground state. This type of conformational relaxation processes coupled with the ligand binding plays is well known for motor proteins motor1 , heme-containing proteins Frauenfelder , etc. In this work, following our general strategy, we want to design a toy ligand with which our toy protein can perform the entire cycle of its ”function.”

We designed a toy protein which has binding site and is able to change its shape. The initial protein conformation and conformational changes activated by the adsorption of the toy ligand are shown in Figure 10. The toy protein was designed as explained above. The linear rearrangement path of the molecule includes $`6`$ conformations. The conformation $`𝒞_1`$ is the lowest energy state for the protein with no ligand. The energies of interaction of the ligand with the ”active center” of the protein are chosen such that as soon as the complex forms, the energies of linear path conformations change in comparison with those in isolated protein molecule. The profile of the energy changes along the path of conformational rearrangements in the protein-ligand complex is shown in Figure 11. When the ligand is attached, the conformation $`𝒞_6`$ of the path has the lowest energy. Hence as the small molecule binds to the protein in the state $`𝒞_1`$, it induces the cascade of conformational changes which drives the protein to the conformation $`𝒞_6`$. The energy profile along the path of the conformational changes corresponds to the case of the suicide inhibitor due to the high barrier for the dissociation of the ligand in the conformation $`𝒞_6`$.

Thus, our model exhibits the conformational response to the ligand binding. When the ligand binds, it jump-starts the cascade of orderly occurring conformational transitions $`𝒞_1𝒞_2\mathrm{}𝒞_6`$. In the computer experiment, we have measured the relaxation times \- the number of Monte Carlo steps in which the toy protein goes from $`𝒞_1`$ into $`𝒞_6`$ upon ligand binding. We repeated this measurement in $`100000`$ independent MC runs, and Figure 12 shows the distribution of relaxation times. The tail of the distribution demonstrates exponential decay that is expected for the biased diffusion in one-dimensional system.

These simulation results are in good agreement with some experimental data. For example, Beece et. al. Frauenfelder  studied the rebinding of carbon monoxide to myoglobin protein after dissociation induced by 1 $`\mu `$s laser impulse. Such a short pulse synchronized conformational transformations of many protein molecules in solution. Under proper circumstances, the population of excited protein molecules exhibited exponential decay, quite similar to our data presented on Figure 12. To reproduce also the multi-exponential or power law decay on the medium time scale observed in some experiments Feng\_Gai  will require some modifications of the present model.

## V Concluding remarks

In this study, we designed a toy model which exhibits quite a few protein-machine-like properties. First and foremost, it has the possibility to move like a mechanical system along the one-dimensional path, which we called rearrangement path. Like real protein, our model can fold into the valley of conformations the bottom of which is the rearrangement path. Thus, it folds, but not only it just folds, it does so to form a functional, movable state. While folded into its low energy valley, our toy protein can move along the path which is the valley bottom. Making only few transverse excursions, it performs a biased random walk along this path. Attachment and detachment of the ligand can switch the direction of the bias, thus making toy protein to go orderly through the closed cycle mimicking protein function.

The usefulness of a toy is the possibility to play. Playing with our toy allows to gain some insights into the general concept of machine-like function of proteins. One question which we can discuss is that of storing energy in the deformed globule. Our model is designed in such a way that the energy of the globule at one end of the rearrangement path is higher than at the other - see, e.g, figure 11, left panel. Can we say this globule acts as a molecular spring? To begin with, this seems to contradict the formulation of the design Hamiltonian which purports to make energy *linear* in the coordinate $`k`$ along the path ($`kW`$), while regular Hookean spring, of course, should have energy quadratic in deformation. In fact, this is unimportant, because design Hamiltonian can be easily modified to purport quadratic (or any other) energy profile; only discrete character of the lattice renders such fine tuning of the model useless. More importantly, all energies involved in our model are in fact free energies. It should be understood that proteins, although machines, are not heat engines: they function in essentially isothermic environment. Therefore, when we say, for instance, that contact of monomers $`i`$ and $`j`$ has energy $`\delta B_{ij}`$, we should have in mind that pre-averaging has been performed over a multitude of small scale, rapidly relaxing degrees of freedom (e.g., $`\chi `$-angles of the residues side groups), and then $`\delta B_{ij}`$ is the free energy of the corresponding contact. From that point of view, we can say that at the upper-energy end of the rearrangement path our model protein stores free energy rather than energy. In other words, it acts not so much like a Hookean spring, but rather like a piece of rubber. This analogy gets even better if we remember that linearity of the strain-stress relation and reversibility of deformation are totally unrelated in the case of rubber. As in other mechanical devices, the question of reversibility in our toy machine is related to its speed: the slower is the rate of operation, the closer it can approach to the limit of being reversible. Of course, when our toy protein undergoes fluctuations along its designed path, it exchanges energy with the surrounding heat bath, because the process occurs at the constant temperature. When conformational relaxation takes place, because random walk along the path is strongly biased - then energy exchange with thermal bath is dominated by the energy transfer into the bath. However, when we imagine that ligand energy has driven the enzyme to the high energy end of its path, then this energy remains available and drives the subsequent conformational relaxation. In this delicate sense, our toy can be said to work like a molecular spring - or, once again, as a piece of rubber. It goes without saying that this is an overdamped spring (or rubber), no oscillations are in question in the course of conformational relaxation.

Our approach is, no doubt, very schematic. It uses to the full extent the discrete geometry of moves on the cubic lattice. This is the heavy price we have to pay for the tractability of the model. In our opinion, the wonderful properties which we found justify the study of this model. We hope also that the fundamental ideas behind our approach will be useful for the more realistic models, including off-lattice ones. These fundamental ideas are the search for a one-dimensional rearrangement pathways in the compact globule and the design of sequences capable of folding into a ”collective funnel” formed together by all states belonging to this pathway. Design of the static protein backbone conformation Tang , as well as sequence design for the static seq\_des\_static  limited flexibility seq\_des\_flex  backbones are all familiar computational approaches. It is a challenge to incorporate the search for the one-dimensional rearrangement pathways into this already difficult area.

## VI Acknowledgements

We acknowledge Minnesota Supercomputing Institute whose computational resources were employed in this study.

AG acknowledges useful discussions and correspondence with L. A. Blumenfeld and A. S. Mikhailov.


format: text

text:
# Renormalization Group Equations in Gauge Theories with Multiple 𝑈(1) Groups

## I Introduction

Renormalization group equations (RGEs) relate parameters at different scales. It plays an important role in the search of physics beyond the standard model (SM) and provides a unique window for physics at extremely high energy. Comprehensive study of two-loop RGEs in a general gauge theory had been performed long ago MV . The calculation was first performed for a simple Lie group. By using suitable substitution rules, these results can be extended to cases of semi-simple Lie groups, as well as cases with only one $`U(1)`$ sector. Recently, these calculations were updated by carefully treating the transformation properties of the fermion fields xiao  and applied to the special case of the SMluo . Extension to cases with multiple $`U(1)`$ groups is more involved, due to the possible presence of kinetic mixing terms between different $`U(1)`$ gauge bosons. In this letter, we extend the calculation to these cases.

Gauge theories with multiple $`U(1)`$ groups can arise in many context. In most grand unified theories, one in general has more than one $`U(1)`$ sector. In string theories, one obtains extra $`U(1)`$ sectors copiously. And there are no general rules to make these extra $`U(1)`$ gauge bosons heavy. Were they discovered, extra $`U(1)`$ bosons would provide important signatures of the underlying theory. In some model buildings, extra gauge $`U(1)`$ bosons were sought to solve the $`\mu `$ problem in supersymmetric theories mu . It is arguably that extra $`U(1)`$ bosons are the most likely new physics next to supersymmetry to be discovered. Over the years, models with extra $`U(1)`$’s have been extensively studied SMU1 .

Unique to multiple $`U(1)`$ gauge theories, there could be kinetic mixing terms in the form

$$\frac{1}{4}\underset{mn}{}\xi _{mn}F^{m\mu \nu }F_{\mu \nu }^n$$
(1)
which are both Lorentz and gauge invariant. One needs extra symmetries to keep them out, while in non-abelian gauge theories, similar terms are forbidden automatically by gauge invariance. Even if one diagonalizes the the abelian gauge fields such that $`\xi _{mn}(mn)`$ vanish at one scale, the mixing terms will in general come back at another scale due to RGE runnings dienes . Especially, these mixings affect the various $`\beta `$\- and $`\gamma `$-functions, which should be included in consistent analysis. Conventionally, one normalizes the abelian gauge fields into canonical forms, such that the kinetic mixing terms vanish at the tree level. For example, the $`\beta `$-functions of gauge couplings have been obtained up to two-loops via this procedure in aguila . However, mixings usually come back in loop-levels and the whole procedure becomes rather complicated. In this letter, we will provide an alternative approach. We retain the kinetic-mixing terms and treat $`\xi _{mn}`$ as free parameters. It proves to be extremely convenient in practice. The $`\beta `$\- and $`\gamma `$-functions of the theory, including the $`\beta `$-functions of the $`\xi _{mn}`$’s, can be obtained in a straightforward manner by extending existing calculations.

The paper is organized as follows, section 2 presents our formalism and the $`\gamma `$-functions of the scalar and fermion fields. Section 3 presents the $`\beta `$-functions of the gauge couplings and those of the $`\xi _{mn}`$’s. Section 4 presents the $`\beta `$-functions of the Yukawa couplings and section 5 those of quartic scalar couplings. We conclude in section 6.

## II formalism

We start with a general renormalizable field theory with a group $`G\times _mU_m(1)`$, where $`G`$ is a compact simple Lie group. The theory contains real scalar fields $`\varphi _a`$ and two-component fermion fields $`\psi _j`$. The Lagrangian can be written as

$$=_0+(\mathrm{gauge}\mathrm{fixing}+\mathrm{ghost}\mathrm{terms})$$
(2)
where

$`_0`$ $`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle \underset{mn}{}}\xi _{mn}F^{m\mu \nu }F_{\mu \nu }^n{\displaystyle \frac{1}{4}}F^{A\mu \nu }F_{\mu \nu }^A`$ (3)
$`+`$ $`{\displaystyle \frac{1}{2}}D^\mu \varphi _aD_\mu \varphi _a+i\psi _j^+\sigma ^\mu D_\mu \psi _j`$
$``$ $`{\displaystyle \frac{1}{2}}(Y_{jk}^a\psi _j\zeta \psi _k\varphi _a+h.c.){\displaystyle \frac{1}{4!}}\lambda _{abcd}\varphi _a\varphi _b\varphi _c\varphi _d,`$
where the diagonal terms with $`\xi _{mm}=1`$ correspond to the usual kinetic terms for abelian gauge bosons, while the off-diagonal $`\xi _{mn}(mn)`$ represent kinetic mixings. The gauge field strengths are

$`F_{\mu \nu }^m=_\mu V_\nu ^m_\nu V_\mu ^m,`$
$`F_{\mu \nu }^A=_\mu V_\nu ^A_\nu V_\mu ^A+gf^{ABC}V_\mu ^BV_\nu ^C,`$
and covariant derivatives of matter fields are

$`D_\mu \varphi _a`$ $`=`$ $`_\mu \varphi _aig\theta _{ab}^AV_\mu ^A\varphi _bi{\displaystyle \underset{m}{}}g_m(Q_m^S)_{ab}V_\mu ^m\varphi _b,`$
$`D_\mu \psi _j`$ $`=`$ $`_\mu \psi _jigt_{jk}^AV_\mu ^A\psi _ki{\displaystyle \underset{m}{}}g_m(Q_m^F)_{jk}V_\mu ^m\psi _k,`$
where $`\theta ^A`$ and $`t^A`$ are the representation matrices of $`G`$ on the scalar and fermion fields, respectively, and $`Q_m^S`$ and $`Q_m^F`$ are their $`U_m(1)`$ charges. Since the scalar fields are real, $`\theta ^A`$ and $`Q_m^S`$ are pure imaginary and antisymmetric. $`Q_m^F`$ are real and diagonal. The constraint on Yukawa coupling matrices imposed by $`U(1)`$ gauge invariance is

$$Y^b(Q_m^S)_{ba}+Y^aQ_m^F+Q_m^FY^a=0,$$
(4)
which can be used to simplify the gauge structures.

Conventionally, one normalizes the $`U(1)`$ gauge fields to assume canonical forms before renormalization. It is easy to show that this normalization does not affect the $`\beta `$-functions of $`g`$, $`Y^a`$ and $`\lambda _{abcd}`$ and the $`\gamma `$-functions of $`V^A`$, $`\varphi _a`$ and $`\psi _j`$ leeyang ; luo2 . On the other hand, the $`\beta `$-functions of $`g_m`$ in the conventional scheme are related to the ones in our scheme by a $`\xi `$-dependent transformation. If one normalizes the $`U(1)`$ gauge fields consistently, both schemes give the same physics prediction luo2 .

In the Landau gauge, the field propagators of the simple gauge group $`G`$ are

$$D_{\mu \nu }^{AB}(k)=\delta ^{AB}\left(g_{\mu \nu }+\frac{k_\mu k_\nu }{k^2}\right)\frac{i}{k^2},$$
(5)
and the propagators of abelian gauge fields are

$`D_{\mu \nu }^{mn}(k)`$ $`=`$ $`\xi _{mn}^1\left(g_{\mu \nu }+{\displaystyle \frac{k_\mu k_\nu }{k^2}}\right){\displaystyle \frac{i}{k^2}},`$ (6)
where $`\xi _{mn}^1`$ is the inverse of $`\xi _{mn}`$. Note that, there are propagations from a $`U_m(1)`$ gauge field to a $`U_n(1)`$ gauge field due to kinetic mixings, with magnitude proportional to $`\xi _{mn}^1`$. For late convenience, we define

$`C_2^{ab}(S)=\theta _{ac}^A\theta _{cb}^A,S_2(S)\delta ^{AB}=\mathrm{Tr}[\theta ^A\theta ^B],`$ (7)
$`C_2^{ab}(F)=t_{ac}^At_{cb}^A,S_2(F)\delta ^{AB}=\mathrm{Tr}[t^At^B],`$ (8)
$`C_2(G)\delta ^{AB}=f^{ACD}f^{BCD},`$ (9)
$`Q_{mn}^{2S}=Q_m^SQ_n^S,Q_{mnpq}^{4S}=Q_m^SQ_n^SQ_p^SQ_q^S,`$ (10)
$`Q_{mn}^{2F}=Q_m^FQ_n^F,Q_{mnpq}^{4F}=Q_m^FQ_n^FQ_p^FQ_q^F,`$ (11)
$`\chi _{mn}=g_mg_n\xi _{mn}^1,`$ (12)
where in Eq. (12), there is no summation over $`m`$ and $`n`$.

Except for the $`\xi _{mn}(mn)`$, we separate the $`\beta `$-function of each coupling constant $`x`$ into two parts,

$$\beta _x=\beta _x(G)+\beta _x(U),$$
(13)
where $`\beta _x(G)`$ is the result when one turns off all $`U(1)`$ interactions and $`\beta _x(U)`$ is the additional contribution when the $`U(1)`$ interactions are turned on. To two-loops, we have

$`\beta _x`$ $`=`$ $`{\displaystyle \frac{1}{16\pi ^2}}\left[\beta _x^{(1)}(G)+\beta _x^{(1)}(U)\right]`$ (14)
$`+`$ $`{\displaystyle \frac{1}{(16\pi ^2)^2}}\left[\beta _x^{(2)}(G)+\beta _x^{(2)}(U)\right].`$
Similarly, we separate the $`\gamma `$-function of field $`i`$

$`\gamma _i`$ $`=`$ $`{\displaystyle \frac{1}{(16\pi ^2)}}\left[\gamma _i^{(1)}(G)+\gamma _i^{(1)}(U)\right]`$ (15)
$`+`$ $`{\displaystyle \frac{1}{(16\pi ^2)^2}}\left[\gamma _i^{(2)}(G)+\gamma _i^{(2)}(U)\right].`$
To two-loops, $`\beta (G)`$ and $`\gamma (G)`$ have been given in MV ; xiao . The additional contributions from $`U(1)`$ groups, $`\beta _x(U)`$ and $`\gamma _i(U)`$ are to be calculated in this letter. These can be obtained from existing calculations by closely inspecting the Feynman diagrams in MV , and making proper substitution of gauge couplings and abelian gauge propagators. The process is straightforward but tedious. We content with by simply reporting the end results here.

### II.1 Scalar wave function renormalization

To one loop, the additional contributions to the $`\gamma `$-functions of scalar fields induced by $`U(1)`$ groups are

$$\gamma _{ab}^{S(1)}(U)=3\chi _{mn}(Q_{mn}^{2S})_{ab},$$
(16)
here and hereafter repeated indices such as $`m`$ and $`n`$ are implicitly summed over all abelian gauge groups, unless specified otherwise. To two loops, the additional contributions are

$`\gamma _{ab}^{S(2)}(U)`$ $`=`$ $`\chi _{mn}\{5\kappa \mathrm{Tr}\left[Q_{mn}^{2F}(Y^{+a}Y^b+Y^{+b}Y^a)\right]`$ (17)
$`+`$ $`{\displaystyle \frac{1}{12}}\chi _{pq}(Q_{mp}^{2S})_{ab}[11\mathrm{T}\mathrm{r}(Q_{nq}^{2S})+40\kappa \mathrm{Tr}(Q_{nq}^{2F})]`$
$`+`$ $`{\displaystyle \frac{3}{2}}\chi _{pq}(Q_{mnpq}^{4S})_{ab}+3g^2\left[C_2(S)Q_{mn}^{2S}\right]_{ab}\},`$
where $`\kappa =\frac{1}{2}(1)`$ for two(four)-component fermions.

### II.2 Fermion wave function renormalization

To one loop, the additional contributions to the $`\gamma `$-functions of fermion fields induced by abelian groups are zero. To two loops,

$`\gamma _F^{(2)}(U)`$ $`=`$ $`\chi _{mn}\{{\displaystyle \frac{1}{4}}Y^aQ_{mn}^{2F}Y^{+a}`$ (18)
$``$ $`{\displaystyle \frac{7}{4}}Q_{mn}^{2F}Y^aY^{+a}+{\displaystyle \frac{9}{2}}(Q_{mn}^{2S})_{ab}Y^aY^{+b}`$
$``$ $`{\displaystyle \frac{1}{4}}\chi _{pq}(Q_{mp}^{2F})[\mathrm{Tr}(Q_{nq}^{2S})+8\kappa \mathrm{Tr}(Q_{nq}^{2F})]`$
$``$ $`{\displaystyle \frac{3}{2}}\chi _{pq}Q_{mnpq}^{4F}.3g^2C_2(F)Q_{mn}^{2F}\}.`$
## III $`\beta `$-functions of gauge couplings and kinetic mixings

### III.1 Nonabelian gauge coupling

At one loop level, the additional contributions to the $`\beta `$-function of gauge coupling $`g`$ of group $`G`$ are zero. To two loop, one has

$$\beta _g^{(2)}(U)=2g^3\chi _{mn}T_{mn}^G(S)+4\kappa g^3\chi _{mn}T_{mn}^G(F)$$
(19)
where $`T_{mn}^G(S)`$ and $`T_{mn}^G(F)`$ are defined as

$`T_{mn}^G(S)\delta ^{AB}`$ $`=`$ $`\mathrm{Tr}(\theta ^A\theta ^BQ_{mn}^{2S})`$
$`T_{mn}^G(F)\delta ^{AB}`$ $`=`$ $`\mathrm{Tr}(t^At^BQ_{mn}^{2F})`$
### III.2 Abelian gauge couplings

The $`\beta `$-function of the gauge coupling $`g_p`$ of group $`U_p(1)`$ is related to the anomalous dimension $`\gamma _p`$ of the corresponding gauge field $`V_\mu ^p`$\[Figure 1(a)\]. The Ward identity ensures that $`\beta _p=g_p\gamma _p`$. To two loops,

$`\beta _p=`$ $``$ $`{\displaystyle \frac{g_p^3}{(4\pi )^2}}\left\{{\displaystyle \frac{4}{3}}\kappa S_2^p(F){\displaystyle \frac{1}{6}}S_2^p(S)+{\displaystyle \frac{2\kappa }{(4\pi )^2}}Y_4^p(F)\right\}`$ (20)
$`+`$ $`{\displaystyle \frac{g_p^3}{(4\pi )^4}}\{2\chi _{mn}[\mathrm{Tr}(Q_{ppmn}^{4S})+2\kappa \mathrm{Tr}(Q_{ppmn}^{4F})]`$
$`+`$ $`2g^2C_2(S)S_2^p(S)+4\kappa g^2C_2(F)S_2^p(F)\}`$
where

$$Y_4^p(F)=\mathrm{Tr}\left[C_2^p(F)Y^aY^{+a}\right],$$
(21)
and there is no summation over $`p`$. $`C_2(R)`$ is the Casimir operator of $`G`$ on the representations $`R`$. $`C_2^p(R)`$ and $`S_2^p(R)`$ are the Casimir operator and the Dykin index of $`U_p(1)`$ acting on fields $`R`$, respectively.

### III.3 Kinetic mixings

The $`\beta `$-functions of the kinetic mixing coefficients $`\xi _{pq}`$ are given by

$$\beta _{pq}=(\gamma _p+\gamma _q)\xi _{pq}+\gamma _{pq}^\xi $$
(22)
where repeated indices don’t imply summation and $`\gamma _{pq}^\xi `$ is the anomalous dimension of the operator $`F^{p\mu \nu }F_{\mu \nu }^q`$\[Figure 1(b)\],

$`\gamma _{pq}^\xi `$ $`=`$ $`{\displaystyle \frac{g_pg_q}{(4\pi )^2}}\left[{\displaystyle \frac{1}{3}}\mathrm{Tr}(Q_{pq}^{2S})+{\displaystyle \frac{8}{3}}\kappa \mathrm{Tr}(Q_{pq}^{2F})\right]`$ (23)
$``$ $`{\displaystyle \frac{g_pg_q}{(4\pi )^4}}\{4\kappa \mathrm{Tr}(Y^aY^{+a}Q_{pq}^{2F}).`$
$`+`$ $`4g^2\left(\mathrm{Tr}\left[Q_{pq}^{2S}C_2(S)\right]+2\kappa \mathrm{Tr}\left[Q_{pq}^{2F}C_2(F)\right]\right)`$
$`+`$ $`4\chi _{mn}(\mathrm{Tr}(Q_{pqmn}^{4S})+.2\kappa \mathrm{Tr}(Q_{pqmn}^{4F}))\}`$
where the indices are summed over $`m(n)`$ but not over $`p(q)`$.

## IV Yukawa coupling

The $`\beta `$-functions of the Yukawa couplings can be expressed as

$$\beta ^a=\gamma ^a+\gamma ^{+F}Y^a+Y^a\gamma ^F+\gamma _{ab}^SY^b$$
(24)
where $`\gamma ^a`$ are the anomalous dimensions of the operators $`\varphi _a\psi _j\zeta \psi _k`$, $`\gamma ^F`$ and $`\gamma _{ab}^S`$ are the anomalous dimensions of the corresponding fermions and bosons, respectively. To one loop level, the additional contributions to the $`\beta `$-function of Yukawa couplings are

$$\beta _a^{(1)}(U)=3\chi _{mn}(Q_{mn}^{2F}Y^a+Y^aQ_{mn}^{2F})$$
(25)
To two loop level, one has

$`\beta _a^{(2)}(U)`$ $`=`$ $`\chi _{mn}\{3\{Q_{mn}^{2F},Y^bY^{+a}Y^b\}.+5Y^b\{Q_{mn}^{2F},Y^{+a}\}Y^b`$ (26)
$``$ $`{\displaystyle \frac{7}{4}}\left[Q_{mn}^{2F}Y^bY^{+b}Y^a+Y^aY^{+b}Y^bQ_{mn}^{2F}\right]`$
$``$ $`{\displaystyle \frac{1}{4}}(Y^bQ_{mn}^{2F}Y^{+b}Y^a+Y^aY^{+b}Q_{mn}^{2F}Y^b)`$
$`+`$ $`6(Q_m^FY^aY^{+b}Q_n^FY^b+Y^bQ_m^FY^{+b}Y^aQ_n^F)`$
$`+`$ $`5\kappa Y^b\mathrm{Tr}\left[Q_{mn}^{2F}(Y^{+a}Y^b+Y^{+b}Y^a)\right]`$
$`+`$ $`6\left[(Q_{mn}^{2S})_{bc}Y^bY^{+a}Y^c2(Q_{mn}^{2S})_{ac}Y^bY^{+c}Y^b\right]`$
$`+`$ $`{\displaystyle \frac{9}{2}}(Q_{mn}^{2S})_{bc}(Y^bY^{+c}Y^a+Y^aY^{+c}Y^b)`$
$``$ $`3g^2\{Q_{mn}^{2F}C_2(F),Y^a\}{\displaystyle \frac{3}{2}}\chi _{pq}\{Q_{mnpq}^{4F},Y^a\}`$
$`+`$ $`6g^2C_2^{ab}(S)\{Q_{mn}^{2F},Y^b\}+6g^2(Q_{mn}^{2S})_{ab}\{C_2(F),Y^b\}`$
$`+`$ $`6\chi _{pq}(Q_{mn}^{2S})_{ab}\{Q_{pq}^{2F},Y^b\}`$
$`+`$ $`{\displaystyle \frac{1}{12}}\chi _{pq}\{Q_{mp}^{2F},Y^a\}[11\mathrm{T}\mathrm{r}(Q_{nq}^{2S})+40\mathrm{T}\mathrm{r}(Q_{nq}^{2F})]`$
$``$ $`21\left[g^2(Q_{mn}^{2S}C_2(S))_{ab}+{\displaystyle \frac{1}{2}}\chi _{pq}(Q_{mnpq}^{4S})_{ab}\right]Y^b`$
$``$ $`{\displaystyle \frac{1}{4}}\chi _{pq}(Q_{mp}^{2S})_{ab}[8\kappa \mathrm{Tr}(Q_{nq}^{2F})+.\mathrm{Tr}(Q_{nq}^{2S})]Y^b\}`$
## V Scalar quartic coupling

The $`\beta `$-functions of the scalar quartic couplings can be expressed as

$$\beta _{abcd}=\gamma _{abcd}+\underset{i}{}\gamma ^S(i)\lambda _{abcd}$$
(27)
where $`\gamma _{abcd}`$ are the anomalous dimensions of the operators $`\varphi _a\varphi _b\varphi _c\varphi _d`$, $`\gamma ^S(i)`$ is the anomalous dimension of the scalar field $`i`$. To one loop level, one has

$`\beta _{abcd}^{(1)}(U)=`$ $``$ $`3\chi _{mn}\mathrm{\Lambda }_{abcd}^{Q,mn}+6g^2\chi _{mn}A_{abcd}^{mn}`$ (28)
$`+`$ $`3\chi _{mn}\chi _{pq}A_{abcd}^{mnpq}`$
and to two loop level, one has

$`\beta _{abcd}^{(2)}(U)`$ $`=`$ $`\chi _{mn}[2\overline{\mathrm{\Lambda }}_{abcd}^{mn}6\mathrm{\Lambda }_{abcd}^{mn}`$
$`+`$ $`4\kappa (H_{abcd}^{Q,mn}H_{abcd}^{F,mn})+5\kappa {\displaystyle \underset{i}{}}Y^{Q2F}(i)\lambda _{abcd}]`$
$`+`$ $`\chi _{mn}\chi _{pq}\{[{\displaystyle \frac{10}{3}}\kappa \mathrm{Tr}(Q_{nq}^{2F})`$
$`+`$ $`{\displaystyle \frac{11}{12}}\mathrm{Tr}(Q_{nq}^{2S})]\mathrm{\Lambda }^{Q,mp}_{abcd}`$
$`+`$ $`{\displaystyle \frac{3}{2}}\mathrm{\Lambda }_{abcd}^{2Q,mnpq}+{\displaystyle \frac{5}{2}}A_{abcd}^{\lambda ,mnpq}+{\displaystyle \frac{1}{2}}\overline{A}_{abcd}^{\lambda ,mnpq}`$
$``$ $`4\kappa (B_{abcd}^{Y,mnpq}10\overline{B}_{abcd}^{Y,mnpq})\}`$
$`+`$ $`g^2\chi _{mn}[3\mathrm{\Lambda }_{abcd}^{SQ,mn}+5A_{abcd}^{\lambda ,mn}`$
$`+`$ $`\overline{A}_{abcd}^{\lambda ,mn}8\kappa (B_{abcd}^{Y,mn}10\overline{B}_{abcd}^{Y,mn})]`$
$``$ $`\chi _{mn}\chi _{pq}\chi _{kl}\{{\displaystyle \frac{15}{2}}A_{abcd}^{mnpqkl}`$
$`+`$ $`[{\displaystyle \frac{7}{3}}\mathrm{Tr}(Q_{np}^{2S})+{\displaystyle \frac{32}{3}}\kappa \mathrm{Tr}(Q_{np}^{2F})]A_{abcd}^{mqkl}\}`$
$``$ $`g^2\chi _{mn}\chi _{pq}\{15A_{abcd}^{Q,mnpq}+{\displaystyle \frac{15}{2}}A_{abcd}^{S,mnpq}`$
$`+`$ $`[{\displaystyle \frac{7}{3}}\mathrm{Tr}(Q_{np}^{2S})+{\displaystyle \frac{32}{3}}\kappa \mathrm{Tr}(Q_{np}^{2F})]A_{abcd}^{mq}\}`$
$``$ $`g^4\chi _{mn}\{{\displaystyle \frac{15}{2}}A_{abcd}^{Q,mn}+15A_{abcd}^{S,mn}`$
$`+`$ $`[{\displaystyle \frac{7}{3}}S_2(S)+{\displaystyle \frac{32}{3}}\kappa S_2(F){\displaystyle \frac{161}{6}}C_2(G)]A_{abcd}^{mn}\}`$
The group factors are defined as

$`\mathrm{\Lambda }_{abcd}^{Q,mn}`$ $`=`$ $`{\displaystyle \underset{i}{}}C_2^{Q,mn}(i)\lambda _{abcd}`$
$`\overline{\mathrm{\Lambda }}_{abcd}^{mn}`$ $`=`$ $`{\displaystyle \frac{1}{8}}{\displaystyle \underset{perms}{}}(Q_{mn}^{2S})_{fg}\lambda _{abef}\lambda _{cdeg}`$
$`\mathrm{\Lambda }_{abcd}^{mn}`$ $`=`$ $`{\displaystyle \frac{1}{8}}{\displaystyle \underset{perms}{}}\lambda _{abef}\lambda _{cdgh}(Q_m^S)_{eg}(Q_n^S)_{fh}`$
$`H_{abcd}^{Q,mn}`$ $`=`$ $`{\displaystyle \underset{i}{}}C_2^{Q,mn}(i)H_{abcd}`$
$`H_{abcd}^{F,mn}`$ $`=`$ $`{\displaystyle \underset{perms}{}}\mathrm{Tr}\left[\{Q_{mn}^{2F},Y^a\}Y^{+b}Y^cY^{+d}\right]`$
$`\mathrm{\Lambda }_{abcd}^{2Q,mnpq}`$ $`=`$ $`{\displaystyle \underset{i}{}}C_2^{Q,mn}(i)C_2^{Q,pq}(i)\lambda _{abcd}`$
$`\mathrm{\Lambda }_{abcd}^{SQ,mn}`$ $`=`$ $`{\displaystyle \underset{i}{}}C_2(i)C_2^{Q,mn}(i)\lambda _{abcd}`$
$`Y^{Q2F}\delta ^{ab}`$ $`=`$ $`\mathrm{Tr}\left[Q_{mn}^{2F}(Y^{+a}Y^b+Y^{+b}Y^a)\right]`$
$`A_{abcd}^{\lambda ,mnpq}`$ $`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle \underset{perms}{}}\lambda _{abef}\{Q_m^S,Q_p^S\}_{ef}\{Q_n^S,Q_q^S\}_{cd}`$
$`A_{abcd}^{\lambda ,mn}`$ $`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle \underset{perms}{}}\lambda _{abef}\{Q_m^S,\theta ^A\}_{ef}\{Q_n^S,\theta ^A\}_{cd}`$
$`\overline{A}_{abcd}^{\lambda ,mnpq}`$ $`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle \underset{perms}{}}\lambda _{abef}\{Q_m^S,Q_p^S\}_{ce}\{Q_n^S,Q_q^S\}_{df}`$
$`\overline{A}_{abcd}^{\lambda ,mn}`$ $`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle \underset{perms}{}}\lambda _{abef}\{Q_m^S,\theta ^A\}_{ce}\{Q_n^S,\theta ^A\}_{df}`$
$`B_{abcd}^{Y,mnpq}`$ $`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle \underset{perms}{}}\{Q_m^S,Q_p^S\}_{ab}\mathrm{Tr}(Q_n^FQ_q^FY^cY^{+d}`$
$`+`$ $`Y^cQ_n^FQ_q^FY^{+d})`$
$`B_{abcd}^{Y,mn}`$ $`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle \underset{perms}{}}\{Q_m^S,\theta ^A\}_{ab}\mathrm{Tr}(Q_n^Ft^AY^cY^{+d}`$
$`+`$ $`Y^cQ_n^Ft^AY^{+d})`$
$`\overline{B}_{abcd}^{Y,mnpq}`$ $`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle \underset{perms}{}}\{Q_m^S,Q_p^S\}_{ab}\mathrm{Tr}(Q_n^FY^cQ_q^FY^{+d})`$
$`\overline{B}_{abcd}^{Y,mn}`$ $`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle \underset{perms}{}}\{Q_m^S,\theta ^A\}_{ab}\mathrm{Tr}(Q_n^FY^ct^AY^{+d})`$
$`A_{abcd}^{mn}`$ $`=`$ $`{\displaystyle \frac{1}{8}}{\displaystyle \underset{perms}{}}\{\theta ^A,Q_m^S\}_{ab}\{\theta ^A,Q_n^S\}_{cd}`$
$`A_{abcd}^{mnpq}`$ $`=`$ $`{\displaystyle \frac{1}{8}}{\displaystyle \underset{perms}{}}\{Q_m^S,Q_p^S\}_{ab}\{Q_n^S,Q_q^S\}_{cd}`$
$`A_{abcd}^{mnpqkl}`$ $`=`$ $`{\displaystyle \underset{i}{}}C_2^{Q,mn}(i)A_{abcd}^{pqkl}`$
$`A_{abcd}^{Q,mn}`$ $`=`$ $`{\displaystyle \underset{i}{}}C_2^{Q,mn}(i)A_{abcd}`$
$`A_{abcd}^{Q,mnpq}`$ $`=`$ $`{\displaystyle \underset{i}{}}C_2^{Q,mn}(i)A_{abcd}^{pq}`$
$`A_{abcd}^{S,mn}`$ $`=`$ $`{\displaystyle \underset{i}{}}C_2(i)A_{abcd}^{mn}`$
$`A_{abcd}^{S,mnpq}`$ $`=`$ $`{\displaystyle \underset{i}{}}C_2(i)A_{abcd}^{mnpq}`$
where $`C_2^{Q,mn}(k)`$ are the eigenvalues of $`Q_m^SQ_n^S`$, and

$`A_{abcd}`$ $`=`$ $`{\displaystyle \frac{1}{8}}{\displaystyle \underset{perms}{}}\{\theta ^A,\theta ^B\}_{ab}\{\theta ^A,\theta ^B\}_{cd}`$
$`H_{abcd}`$ $`=`$ $`{\displaystyle \underset{perms}{}}{\displaystyle \frac{1}{4}}\mathrm{Tr}(Y^aY^{+b}Y^cY^{+d})`$
## VI Conclusion

We have now presented the two-loop RGEs for all dimensionless parameters in a general gauge theory with multiple abelian gauge groups. We have retained all kinetic mixing terms and did not canonically normalize the abelian gauge fields. This approach proved to be extremely convenient. Were all abelian gauge interactions turned off, the two-loop RGEs were given in MV ; xiao . When the abelian gauge interactions are turned on, the additional contributions have been obtained from existing calculations by a close inspection of the relevant Feynman diagrams and by using suitable substitution rules. Were $`G`$ not simple, but rather a direct product of simple groups $`G_k`$ with gauge coupling $`g_k`$, the RGEs can again be obtained by using substitution rules as given in MV ; xiao . For part of $`g^4`$-related terms, the following extra substitution rules are needed,

$`g^4A_{abcd}^{Q,mn}{\displaystyle \frac{1}{8}}g_k^2g_l^2{\displaystyle \underset{i}{}}C_2^{Q,mn}(i){\displaystyle \underset{perms}{}}\{\theta _k^A,\theta _l^B\}_{ab}\{\theta _k^A,\theta _l^B\}_{cd}`$
$`g^4A_{abcd}^{S,mn}{\displaystyle \frac{1}{8}}g_k^2g_l^2{\displaystyle \underset{i}{}}C_2^k(i){\displaystyle \underset{perms}{}}\{\theta _l^A,Q_m^S\}_{ab}\{\theta _l^A,Q_n^S\}_{cd}`$
$`g^4S_2(R)A_{abcd}^{mn}{\displaystyle \frac{1}{8}}g_k^4S_2^k(R){\displaystyle \underset{perms}{}}\{\theta _k^A,Q_m^S\}_{ab}\{\theta _k^A,Q_n^S\}_{cd}`$
$`g^4C_2(G)A_{abcd}^{mn}{\displaystyle \frac{1}{8}}g_k^4C_2^k(G){\displaystyle \underset{perms}{}}\{\theta _k^A,Q_m^S\}_{ab}\{\theta _k^A,Q_n^S\}_{cd}`$
By introducing a dummy field, $`\beta `$-functions of the dimensional parameters can also be deduced from those of $`Y^a`$ and $`\lambda _{abcd}`$ xiao ; luo2 .

###### Acknowledgements.

The work was supported by a Fund for Trans-Century Talents, CNSF-90103009 and CNSF-10047005.