Thesis (Index)   <-  Sean Forman   <-  You Are Here

Tweaking Algorithms for $\beta$-Strand Alignment

In the course of folding a protein with a discrete set of backbone torsion angles, HOPS often creates $\beta$-strands which do not take part in hydrogen bonding with other nearby $\beta$-strands to form $\beta$-sheets. These unformed $\beta$-sheets are poorly scoring conformations. However, with a small change or ``tweak'' to a handful of amino acids, the $\beta$-strands could be brought into alignment and the resulting conformation would have a significantly better score.

Discretization of a protein's $\phi - \psi$ torsion angles and the generation of a finite number of folds converts an overwhelming continuous problem into a form slightly easier, but still NP-hard and quite difficult to solve. This simplification, however, must come with trade-offs. What we find using our technique is that isolated ``local'' motifs occur with great regularity: $\alpha$-helices and $\beta$-strands. However, these sheet structures do not routinely come together to form more common super-secondary structures like $\beta$-sheets, or any sort of aligned helical structures. These structures are energetically advantageous, but they involve complicated long-range interactions which are unlikely to occur in a discrete search space. This isn't just a shortcoming of our technique; most ab initio techniques perform much better on helical proteins than on those with mostly $\beta$-sheets.

This stands to reason. These structures are largely global in nature. The alignment of the beginning of a $\beta$-strand, A, to the end of $\beta$-strand, B, can be achieved in a vast number of ways. In nature, you can think of the two strands arriving close together and then pulling each other into alignment. When pulled into alignment, the turn regions between the two $\beta$-strands are changed slightly to bring them together, so hydrogen bonding can occur. This is a relaxation of the discrete $\phi - \psi$ angle pairs that we are using.

Due to this, the alignment of the $\beta$-strands defies discretization. Certainly it would be possible to identify, describe, and implement some canonical set of turn angles. This would aid in the creation of anti-parallel $\beta$-sheets, but this is in some sense is the least interesting behavior. A far more interesting case is where two strands from distant regions of the protein come together to form either a parallel or anti-parallel $\beta$-sheet (see Figure 5.1).

[ Sheet structure for protein G.] Two ``local'' anti-parallel $\beta$-sheets come together to form a parallel sheet. Amino acid eight and amino acid fifty-six are taking part in a hydrogen bond. An example of global behavior

In order to capture this global nature of secondary structure, we need to shift from discrete torsion angles to continuous torsion angles. While the discrete angles may bring two structures near each other, they will not close the last few angstroms of gap and properly align the two protein segments. Ideally, we would take the discrete angles and alter them slightly in a continuous manner in order to align two $\beta$-strands. This means that we will not be making large changes in the protein structure-only small tweaks.

This chapter first gives sufficient background to understand the geometric and mathematical techniques HOPS uses to build protein conformations. This is followed by a description of an optimization technique which relaxes the turn region between two $\beta$-strands in order to create a $\beta$-sheet.

Homogeneous Coordinate Systems

Our representation of the protein is based on homogeneous coordinate systems [19]. In robotics, homogeneous coordinate systems are used to model the position of a robot arm. The robot arm (or more generally a rigid link open kinematic chain) typically has several rigid pieces connected by linkages. These linkages have some number of degrees of freedom that allow them to change the arm's configuration. Beginning with the position of the base, one can progressively compute the position of the robot's effector (the business end of the arm) by computing the positions of the intermediate pieces using the position of the linkages and the lengths of the arms.

Similarly, our protein backbone is a set of fixed bond lengths and fixed bond angles but variable torsion angles. Using the same idea as with robotics, we are able to find the global position of all the backbone atoms beginning with the position of the first amino acid's $\mathbf{N}$ atom through to the final amino acid's $\mathbf{C'}$ atom (see Figure 5.2).

Figure 5.2: Protein backbones and robot arms can be modeled in very similar ways.

Three-dimensional geometry has six degrees of freedom: the atom's 3-D location, $(x,y,z)$, and also the atom's orientation (its yaw, pitch and roll, or more precisely the angle of rotation around the $z, y$, and $x$ axes). To simplify the manipulation of these positions and the calculation of subsequent atomic positions, these six degrees of freedom are expanded into twelve entries in a $4 \times 4$ matrix. While this takes twice as much storage space, it reduces the complexity of updating the atomic positions substantially.

We need two types of $4 \times 4$ matrices in order to compute the position of each backbone atom in this manner. Each backbone atom has its own global position matrix, and each covalent bond connecting two backbone atoms has a related transform matrix. As we travel from backbone atom A to adjoining backbone atom B, we compute atom B's position matrix by multiplying atom A's position matrix by the transform matrix of the bond connecting them. Using this technique, we can find the global location of every atom in the protein using only the torsion angle and the angle and length of the bonds between the backbone atoms along with the position of the first backbone atom. This first backbone atom can be set to any arbitrary global position,⁴¹but we typically set the first backbone atom at the origin and the first backbone bond oriented along the global $x$-axis.

Equation 5.1 shows how we can convert $\left(l_X, l_Y, l_Z\right)$, the three location coordinates, and the rotations about the three axes, $\left(r_X, r_Y, r_Z\right)$, into a position matrix. The order in which the rotations are preformed is important, and in the equation given, the first rotation is about the $x$-axis, followed by $z$ and $y$. The atom's location can always be read in the first three rows of the fourth column.

$$\begin{multline} \mathnormal{P} \left(r_X, r_Y, r_Z, l_X, l_Y, l_Z \right) = ... ...cos r_X \cos r_Y & l_Z \\ 0 & 0 & 0 & 1 \\ \end{array}\right) \end{multline}$$

As we build the protein, this equation is not actually used other than when initializing the very first backbone atom. This position matrix is set to the $\mathnormal{P}_{1,\ensuremath{\mathbf{N}}} (0,0,0,0,0,0)$. Note that our subscript for $\mathnormal{P}$ tells us that this is the position matrix for the first amino acid's $\mathbf{N}$ atom.

$$\mathnormal{P}_{1,\ensuremath{\mathbf{N}}} = \left( \begin{array}... ...0 0 & 1 & 0 & 0 0 & 0 & 1 & 0 0 & 0 & 0 & 1 \end{array} \right)$$ (12)

As stated before, following backbone atom position matrices will be the product of this initial position matrix and all intervening transform matrices. Just as with the position matrices the transform matrix is a $4 \times 4$ matrix. In fact, there is a transform matrix for each of the six basic transformations we can make to an atom's position. We can rotate the atom's position about any of the three axes, or we can translate it along any of its three axes. If we have a series of these basic transforms to perform, we can do that by multiplying the basic matrices in the order in which the transformations occur. In our case, we translate the length of the bond along the $x$-axis, followed by a rotation about the $x$-axis, which is the torsion angle for that bond, and then rotate about the $z$-axis, which is the bond angle. Here $\mathbf{T}$ is the transform matrix that represents a bond with a length, $L$, a bond angle of $\kappa$ and a torsion angle of $\tau$ (see Equation 5.3). Just as we have position matrices for each angle, we can define transform matrices for each bond, $\mathnormal{T}_{i,\theta}$, where this is the transform matrix for amino acid $i$'s $\theta$ bond ( $\theta \in \{\phi,\psi,\omega\}$). This converts the bond's three degrees of freedom to nine entries in a $4 \times 4$ matrix.

$$\mathnormal{T} \mathbf{ \left(\tau, \kappa, \mathnormal{L} \right... ...\cos \kappa \sin \tau & \cos \tau & 0 \\ 0 & 0 & 0 & 1 \\ \end{array} \right)$$ (13)

As stated, the computation of a backbone atom's position matrix requires the previous atom's position matrix and the transformation matrix representing the bond between them. The position matrices for amino acid $i$ can be found if the position matrix for $\ensuremath{\mathbf{C^\prime}}$ of amino acid $i-1$ is known and the three transformation matrices corresponding to the $\omega,~\phi,$ and $\psi$ torsion angles can be found. Figure 5.3 shows how the coordinate system changes as we progress along the protein backbone, and equations 5.4 through 5.6 show how we compute the positions of the backbone atoms as we move along the protein.

$$\mathnormal{P}_{i,\ensuremath{\mathbf{N}}} = \mathnormal{P}_{i - 1,\ensuremath{\mathbf{C'}}} \mathnormal{T}_{i - 1,\omega}$$ (14)

$$\mathnormal{P}_{i,\ensuremath{\mathbf{C_{\alpha}}}} = \mathnormal{P}_{i,\ensuremath{\mathbf{N}}} \mathnormal{T}_{i,\phi}$$ (15)

$$\mathnormal{P}_{i,\ensuremath{\mathbf{C'}}} = \mathnormal{P}_{i,\ensuremath{\mathbf{C_{\alpha}}}} \mathnormal{T}_{i,\psi}$$ (16)

$$\mathnormal{P}_{i,\ensuremath{\mathbf{C'}}} = \mathnormal{P}_{i -... ...} \mathnormal{T}_{i - 1,\omega} \mathnormal{T}_{i,\phi} \mathnormal{T}_{i,\psi}$$ (17)

[A short backbone snip with the coordinate systems representing potential orientation matrices overlaid.] The orientation changes as we move along the backbone. In addition to the rotation of the axes, the origin of the coordinate system translates from position $A$ to position $B$

As Equation 5.7 shows, the effect of the transform matrices is cumulative. Knowing only the bond angle and length and the torsion angle of each bond, makes it possible for us to find the position of the final atom in the protein backbone. Later on when we align two $\beta$-strands, we will be altering the $\mathnormal{T}_{i,\theta}$ matrices slightly in order to change the position of a set of atoms.

Tweaking

In Equation 5.3, we saw that the transform matrix from one bond to another is determined by the bond angle and length and the bond's torsion angle. Since the first two parameters are essentially constant, we find that the transform matrices are matrices of one variable, the torsion angle for the bond. This also points out that altering a single torsion angle propogates the change throughout the entire protein.⁴²

This means that, at least in our algorithm, the set of torsion angles chosen to build the protein uniquely define its structure. In the case of real proteins, the bond lengths and angles do vary, but not to a great degree. As stated earlier (see Section 2.1.1), these variations can be determined using experimental methods and are known to be small.

Formulation of the Tweaking Problem

Let's for a moment assume we have a $\beta$-strand A and another $\beta$-strand B with a set of amino acids $\{a_1, ..., a_n \}$ between them. Using our discrete angles we may have something like Figure 5.4.

Figure 5.4: A pair of misaligned $\beta$-strands.

A and B are too far apart to bond with each other as the current configuration stands. This is likely to be a poor partial fold using any common energy function.⁴³ In HOPS, such a partial fold when found would be pruned. However, some small changes to the backbone torsion angles can lead to a conformation like Figure 5.5. This conformation will be much better energetically and would be more likely to be part of a final prediction.

Figure 5.5: A pair of aligned $\beta$-strands.

This points out that the difference between a bad conformation and a good one is often small; the ability to ``relax'' the backbone torsion angles of a turn region can cause sizeable improvements in the energy of a protein. To do this, we need a way to align a pair of amino acids and, then be certain that a $\beta$-sheet will then propogate as new amino acids are added.

Structures like $\alpha$-helices and $\beta$-strands are largely local in nature. Given the proper pair of $\phi - \psi$ torsion angles to go along with an $\omega$ of $180^{\circ}$ we can set up the case where the $\alpha$-helix or $\beta$-strand propogates infinitely. Using this fact for $\beta$-strands, we can define the position of A's $\mathbf{C'}$ relative to B's $\mathbf{N}$ in some way such that the two strands will align no matter how long of a $\beta$-sheet we decide to build. (It may be helpful to review the bonding pattern of $\beta$-sheets from Section 2.1.3.)

Given this ability, our problem is no longer to align two strands, but now to align two atoms.

If you recall in Section 5.1, we described the manner in which the global position of each amino acid was determined. More specifically, the position of our backbone atoms can be found by multiplying the $4 \times 4$ position matrix, $\mathnormal{P}$, of the previous backbone atom with the $4 \times 4$ transform matrix, $\mathnormal{T}$, representing the bond between the two atoms. $\mathnormal{T}$ depends solely upon the bond length, the bond angle and torsion angle. The bond angle and length are constants, and we only vary the torsion angles, $\tau_{i,\theta}$, so this becomes a matrix of only one variable (see Equation 5.3).

Since we only need to set up one interaction between the two $\beta$-strands, we can simplify the problem. We only look at aligning amino acid $i$'s $\mathbf{C'}$ with amino acid $j$'s $\mathbf{N}$. The tweakable torsion angles are then assumed to be the backbone torsion angles for the amino acids $i+1$ to $j-1$. However, we can create cases where only some subset of the intervening amino acids will be tweakable. This might occur if a $\alpha$-helix occurs between the two $\beta$-strands we are attempting to align. In this case, we don't want to alter the torsion angles within the $\alpha$-helix as this may destroy the helical structure. These intervening structures are considered rigid bodies unchangeable by our tweaking algorithm.

Our problem is then: given a position matrix $\mathnormal{P}_{j,\ensuremath{\mathbf{C^\prime}}}$ and a desired (or goal) position matrix of $\mathnormal{P}_{k,\ensuremath{\mathbf{N}}}$ find a set of transform matrices, $\{\mathnormal{T}_{i,\theta}\}, i = j+1, ..., k - 1, \theta = \{\phi,\psi,\omega\}$ such that

$$\mathnormal{P}_{j,\ensuremath{\mathbf{C^\prime}}} \mathnormal{T}_... ..... \mathnormal{T}_{k-1,\psi} \mathnormal{T}_{k-1,\omega} = \mathnormal{P}_{k,N}$$ (18)

Each of these transform matrices is a function of the torsion angle and any changes made to it. We signify this as, $${ \mathnormal{T}_{i,\theta} ( \tau_{i,\theta} + \Delta \tau_{i,\theta})}$$. This means that if we find $\{ \Delta \tau_{i,\theta} \}$ that equate the left and right hand sides of Equation 5.8 then the two atoms ( $\ensuremath{\mathbf{C'}}_j$ and $\ensuremath{\mathbf{N}}_k$) will be in proper alignment for the creation of a $\beta$-sheet.

As mentioned earlier (Section 5.1) each $4 \times 4$ position matrix is an expansion of only six variables (three location and three orientation). Expecting each entry of the two matrices in Equation 5.8 to match exactly would overdefine the problem. Instead, we can solve this problem by requiring the two matrices to agree on just six entries, the rightmost column of locations, $[l_x l_y l_z 1]^T$, and three entries from the upper left corner $3 \times 3$ sub-matrix. We found that using the $(1,1)$, $(2,1)$, and $(3,3)$ entries was sufficient to force the three orientation values to be equal. In other words, if we know the indicated entries from the atom's $4 \times 4$ position matrix we know the atom's position,

$$\left( \begin{matrix}\times & - & - & \times \times & - & - & \times - & - & \times & \times - & - & - & - \end{matrix} \right)$$ (19)

and if two matrices agree on those six entries, they describe the same atom position. We will use the symbol, $\doteqdot$, to represent this relation. This is a weaker relation than matrix equality, so all the properties true in matrix equality hold here as well.

We can further simplify the equation by removing $\mathnormal{P}_{j,\ensuremath{\mathbf{C^\prime}}}$ and $\mathnormal{P}_{k,\ensuremath{\mathbf{N}}}$ from the picture. If we left-multiply both sides by $\mathnormal{P}_{j,\ensuremath{\mathbf{C^\prime}}}^{-1}$, we are left with

$$\mathnormal{T}_{j,\omega} \mathnormal{T}_{j+1,\phi} \mathnormal{T... ...^{-1} \mathnormal{P}_{k,\ensuremath{\mathbf{N}}} \doteqdot \mathnormal{G_\zeta}$$ (20)

where $\zeta$ is either parallel ($para$) or anti-parallel ($anti$), and $\mathnormal{G_\zeta}$ is a canonical position matrix for the creation of a $\beta$-sheet. The values for $\mathnormal{G_\zeta}$ are derived from observations of actual $\beta$-sheets and experiments with efforts to align $\beta$-sheets in HOPS.

$$\mathnormal{G}_{anti} = \left( \begin{matrix}-0.991 & - & - & -2.... ... - & - & 3.250 - & - & 1.000 & 0.000 - & - & - & - \end{matrix} \right)$$ (21)

$$\mathnormal{G}_{para} = \left( \begin{matrix}0.500 & - & - & -1.4... ... & - & 3.930 - & - & -1.000 & -0.140 - & - & - & - \end{matrix} \right)$$ (22)

We have one more constraint to make. We want as much continuity as possible from the initial set of torsion angles to the new tweaked set. A wildly different set of torsion angles is more likely to lead to a steric clash within the turn region, and to be true to our discretization we wish to choose angles that are relatively close to the discrete angles that are the basis of our search space.

Our problem can now be formulated as a least-squares equality constrained problem.

Solution of the Tweaking Problem

Defined formally, we want to find a solution to the following problem.

$$\operatorname{minimize} \qquad \sum_{(i,\theta) = (j+1,\phi) }^{ (k-1,\omega)} \Vert \Delta \tau_{i,\theta} \Vert _2$$ (23)

$$\operatorname{subject to} \qquad \mathnormal{T}_{j,\omega} ( \tau_{j,\omega}) \prod_{(i,\th... ...eta} ( \tau_{i,\theta} + \Delta \tau_{i,\theta}) \doteqdot \mathnormal{G_\zeta}$$

or restated using vector notation as

$$\operatorname{minimize} \qquad \Vert \boldsymbol{\Delta \tau} \Vert _2$$ (24)

$$\operatorname{subject to} \qquad \mathbf{T} ( \boldsymbol{\tau + \Delta \tau}) \doteqdot \mathnormal{G_\zeta}$$ (25)

where $\mathbf{T}$ is the product of the transform matrices for the amino acids between the two that we are attempting to align, and $\boldsymbol{\tau}$ and $\boldsymbol{\Delta \tau}$ are vectors representing the original torsion angles and the amount of change in the torsion angles, respectively. We are holding the torsion angle for the $${\mathnormal{T_{j,\omega}}}$$ transform matrix constant. This is done so that we do not unnecessarily increase the chance of a steric clash between the amino acids we are attempting to align. We could invert this matrix and multiply it with $\mathnormal{G_\zeta}$. Since we are going to convert the main constraint into a linear function anyway, it will only affect the coefficient matrix for our system of linear equations, so we will fold it into the main constraint function.

This problem as stated is difficult to solve optimally, as $\mathbf{T}$ is non-linear. Therefore, we must settle for an approximate solution which we improve iteratively. While we will make certain any solution satisfies our equality constraint, we will not be certain that the objective function is fully minimized. For our purposes, this is sufficient as our overarching goal is to satisfy the equality constraint which will align the two $\beta$-strands.

To solve this problem, we need to convert the constraint into a linear function. This is done easily using a Taylor polynomial approximation of the main constraint.

$$\mathbf{T} (\boldsymbol{\tau + \Delta \tau}) \thickapprox \mathbf... ...d\mathbf{T}(\boldsymbol{\tau_{}})}{d\boldsymbol{\tau}}}\boldsymbol{\Delta \tau}$$ (26)

The use of an approximation does not guarantee that we satisfy Equation 5.15. And we will actually have two severely mis-aligned $\beta$-strands if we were to just solve this system of equations once and then apply our solution to the protein. Instead, we iteratively update the quality of our Taylor approximation with solution from the previous iteration. As the quality of this approximation improves, we get closer to the ultimate goal of $${\mathbf{T} ( \boldsymbol{\tau + \Delta \tau}) \doteqdot \mathnormal{G_\zeta}}$$.

To incorporate this iterative technique, we define $${ \boldsymbol{\tau_{i + 1}} = \boldsymbol{\tau_0 + \Delta \tau_{i + 1}}}$$. This is slightly different than techniques where each solution is a function of the previous iteration's solution. In our case, we use the previous iteration to set up the problem to solve. Our problem now becomes.

$$\operatorname{minimize} \qquad f(\boldsymbol{\Delta \tau_{i+1}}) = \tfrac{1}{2} \Vert \boldsymbol{\Delta \tau_{i+1}} \Vert _2^2$$ (27)

$$\operatorname{subject to} \qquad \mathbf{T} ( \boldsymbol{\tau_{i+1}}) \doteqdot \mathnormal{G_\zeta}$$ (28)

We have slightly changed the objective function $f$ so that $\triangledown f = \sum \Delta \tau$.⁴⁴Next, we perform a series of operations on (5.18) in order to satisfy it iteratively.

$$\mathbf{T} ( \boldsymbol{\tau_{i+1}}) \doteqdot \mathnormal{G_\zeta}$$ (5.18)

$$\mathbf{T} ( \boldsymbol{\tau_{i} + \tau_{i+1} - \tau_{i}}) \doteqdot \mathnormal{G_\zeta}$$

$$\mathbf{T} ( \boldsymbol{\tau_{i} + (\Delta \tau_{i+1} - \Delta \tau_{i})}) \doteqdot \mathnormal{G_\zeta} \quad(\boldsymbol{\tau_{i + 1}} = \boldsymbol{\tau_0 + \Delta \tau_{i + 1}})$$

$$\mathbf{T} ( \boldsymbol{\tau_{i}}) + \ensuremath{\frac{d\mathbf{... ...ldsymbol{\tau}}}(\boldsymbol{\Delta \tau_{i+1}} - \boldsymbol{\Delta \tau_{i}}) \doteqdot \mathnormal{G_\zeta} \quad( )$$

Solving for $\boldsymbol{\Delta \tau_{i+1}}$, we now have a linear constraint in terms of $\boldsymbol{\Delta \tau_{i+1}}$. In the previous iteration, we found $\boldsymbol{\Delta \tau_{i}}$ and $\boldsymbol{\tau_{i}}$.

$$\ensuremath{\frac{d\mathbf{T}(\boldsymbol{\tau_{i}})}{d\boldsymbo... ...hbf{T}(\boldsymbol{\tau_{i}})}{d\boldsymbol{\tau}}}\boldsymbol{\Delta \tau_{i}}$$ (29)

and our least-squares linear equality-constrained program now becomes

$$\operatorname{minimize} \qquad f(\boldsymbol{\Delta \tau_{i+1}}) = \tfrac{1}{2} \Vert \boldsymbol{\Delta \tau_{i+1}} \Vert _2^2$$ (30)

$$\operatorname{subject to} \qquad \ensuremath{\frac{d\mathbf{T}(\boldsymbol{\tau_{i}})}{d\bo... ...hbf{T}(\boldsymbol{\tau_{i}})}{d\boldsymbol{\tau}}}\boldsymbol{\Delta \tau_{i}}$$ (31)

where $${\boldsymbol{\Delta \tau_{i+1}} = \boldsymbol{\tau_{0}} + \boldsymbol{\Delta \tau_{i+1}}}$$ and $${\boldsymbol{\Delta \tau_{0}} = \mathbf{0}}$$. We continue iterating until

$$\displaystyle{ \Vert \mathbf{T} (\boldsymbol{\tau_{i+1}}) - \mathnormal{G_\zeta} \Vert < \epsilon}$$ (32)

or it becomes apparent that we will not be able to find a solution. See Section 5.2.3 for a more detailed discussion of termination criteria.

In our main constraint (5.21), we will see that $${\ensuremath{\frac{d\mathbf{T}(\boldsymbol{\tau_{}})}{d\boldsymbol{\tau}}}}$$ is built as a $6 \times m$ matrix, where $m$ is the number of tweakable backbone torsion angles between the two amino acids we wish to align, $\boldsymbol{\Delta \tau}$ is an $m \times 1$ array, and the right-hand side can be written as a $6 \times 1$ array.⁴⁵ This will allow us to replace $\doteqdot$ with an equal sign ($=$).

Prior to building this equality constraint, we will need to define a transform derivative matrix ( $${d\mathnormal{T \left(\tau \right) }}$$). If you recall the matrix introduced in Equation 5.3 which computes our transformation matrix, we see that the matrix is a function of $\tau_{i,\theta},~\kappa_{i,\theta},$ and $L$. However, we treat the bond angle, $\kappa$, and the bond length, $L$ as constants for every backbone bond, so this matrix has just one variable the torsion angle, $\tau$. If we take the derivative of the matrix with respect to our lone variable we get the derivative matrix, $d\mathnormal{T_{i,\theta}}$.

$$d\mathnormal{T_{i,\theta}} = \left( \begin{array}{rrrr} 0 & 0 & 0... ...cos \kappa \cos \tau & -\sin \tau & 0 0 & 0 & 0 & 0 \end{array} \right)$$ (33)

With this machinery we can then build both sides of our equality constraint from Equation 5.33. In the definition of $${\ensuremath{\frac{d\mathbf{T}(\boldsymbol{\tau_{}})}{d\boldsymbol{\tau}}}}$$, each column corresponds to the partial derivative of $\mathbf{T}$ with respect to a specific $\tau_{i,\theta}$. Since $\mathbf{T}$ is a product of many transform matrices, $${\mathnormal{T}_{i,\theta}}$$, each of which is a function of a single variable, we can simply apply the product rule to compute the partial derivative.

$$\frac{\partial \mathbf{T}}{\partial \tau_{i,\theta}} = \mathnorma... ...hnormal{T}_{i,\theta} ... \mathnormal{T}_{k-1,\psi} \mathnormal{T}_{k-1,\omega}$$ (34)

Then we build $${\ensuremath{\frac{d\mathbf{T}(\boldsymbol{\tau_{}})}{d\boldsymbol{\tau}}}}$$ by computing $${\frac{\partial \mathbf{T}}{\partial \tau_{i,\theta}}}$$ for each backbone bond we are allowing to be tweaked, and then using the six standard entries (5.9) from our $\doteqdot$ relationship as a single six-entry column. Applying this for each of the $m$ $\tau_{i,\theta}$ values gives us the full matrix and defines everything in our least-squares equality-constrained problem (5.21).

In order to solve each instance of our program, we apply Lagrangian multipliers to the problem, which produces a system of linear equations whose solution will solve the given program (5.20).

$$\mathbf{I}\boldsymbol{\Delta \tau}_{i + 1} \quad- \quad\ensuremath{\frac{d\mathbf{T}(\boldsymbol{\tau_{i}})}{d\boldsymbol{\tau}}}^\mathnormal{T} \boldsymbol{\lambda} \quad= \quad\mathbf{0}$$ (35)

$$\ensuremath{\frac{d\mathbf{T}(\boldsymbol{\tau_{i}})}{d\boldsymbol{\tau}}}\boldsymbol{\Delta \tau}_{i + 1} \quad- \quad\boldsymbol{0 \lambda} \quad= \quad\mathnormal{G_\zeta} - \mathbf{T} (\boldsymbol{\tau}_{i}) + ... ...hbf{T}(\boldsymbol{\tau_{i}})}{d\boldsymbol{\tau}}}\boldsymbol{\Delta \tau_{i}}$$

If we let $\mathbf{A} = \displaystyle{\ensuremath{\frac{d\mathbf{T}(\boldsymbol{\tau_{i}})}{d\boldsymbol{\tau}}}}$ and $${\mathbf{b} =\mathnormal{G_\zeta} - \mathbf{T} (\boldsymbol{\tau}_{i}) + \mathbf{A} \boldsymbol{\Delta \tau_{i}} }$$, then we have the following system.

$$\mathbf{I}\boldsymbol{\Delta \tau}_{i + 1} \quad- \quad\mathbf{A}^\mathnormal{T} \boldsymbol{\lambda} \quad= \quad\mathbf{0}$$ (36)

$$\mathbf{A} \boldsymbol{\Delta \tau}_{i + 1} \quad- \quad\boldsymbol{0 \lambda} \quad= \quad\mathbf{b}$$ (37)

Solving this system utilizes some standard techniques from numerical linear algebra. For reasons that will become apparent we are going to operate on $\mathbf{A}^\mathnormal{T}$, which is $m \times 6$ and is a real matrix. For any real matrix, we are able to produce a factorization, $\mathbf{A}^\mathnormal{T} = \mathbf{QR}$ where $\mathbf{Q}$, which is $m \times 6$, is an orthogonal matrix and $\mathbf{R}$, which is $6 \times 6$, is upper triangular. A matrix is orthogonal if $\mathbf{QQ}^\mathnormal{T} = \mathbf{I}$ and $\mathbf{Q}^\mathnormal{T}\mathbf{Q} = \mathbf{I}$, and a matrix is upper triangular if it is square and has only zero entries below the diagonal. A system of linear equations with an upper triangular matrix, $\mathbf{Rx} = \mathbf{b}$ can be solved much more quickly than a full matrix can be solved [3]. This difference is not insignificant. Using Gaussian Elimination on the full matrix, we would take many more operations to form an upper triangular matrix ($O(n^3)$) and then solve the system of equations. The factorization allows us to concentrate on an $m \times 6$ matrix rather than a $(m+6) \times (m+6)$ matrix. Using our factorization, $\mathbf{A}^\mathnormal{T} = \mathbf{QR}$ and $\mathbf{A} = \mathbf{R}^\mathnormal{T}\mathbf{Q}^\mathnormal{T}$ Equation (5.26) becomes,

$$\mathbf{I}\boldsymbol{\Delta \tau}_{i + 1} \quad \quad- \quad\mathbf{QR} \boldsymbol{\lambda} \quad \quad= \quad\mathbf{0}$$ (38)

$$\mathbf{A}\mathbf{I}\boldsymbol{\Delta \tau}_{i + 1} \quad \quad- \quad\mathbf{A}\mathbf{QR} \boldsymbol{\lambda} \quad \quad= \quad\mathbf{A}\mathbf{0}$$ (39)

$$\mathbf{A} \boldsymbol{\Delta \tau}_{i + 1} \quad \quad- \quad\mathbf{R}^\mathnormal{T}\mathbf{Q}^\mathnormal{T}\mathbf{QR} \boldsymbol{\lambda} \quad \quad= \quad\mathbf{0}$$ (40)

$$\mathbf{A} \boldsymbol{\Delta \tau}_{i + 1} \quad \quad- \quad\mathbf{R}^\mathnormal{T}\mathbf{R} \boldsymbol{\lambda} \quad \quad= \quad\mathbf{0}$$ (41)

$$\mathbf{A} \boldsymbol{\Delta \tau}_{i + 1} \quad \quad \quad= \quad\mathbf{R}^\mathnormal{T}\mathbf{R} \boldsymbol{\lambda}$$ (42)

Then substituting this into (5.27), we have a simplified system to solve.

$$\quad\mathbf{R}^\mathnormal{T}\mathbf{R} \boldsymbol{\lambda} = \mathbf{b}$$ (43)

If we let $${\mathbf{R} \boldsymbol{\lambda} = \mathbf{z}}$$, we now have a sequence of three systems of linear equations that we can use to solve for $${\boldsymbol{\Delta \tau}_{i + 1}}$$.

$$\mathbf{R}^\mathnormal{T} \mathbf{z} = \mathbf{b}$$ (44)

$$\mathbf{R} \boldsymbol{\lambda} = \mathbf{z}$$

$$\boldsymbol{\Delta \tau}_{i + 1} = \mathbf{A}^\mathnormal{T} \boldsymbol{\lambda}$$

The three systems are simple to solve because $\mathbf{R}$ is upper triangular (and $\mathbf{R}^\mathnormal{T}$ is lower triangular) which makes $\boldsymbol{\lambda}$ and $\mathbf{z}$ easy to compute, and then $\boldsymbol{\Delta \tau}_{i + 1}$ is computed using just a simple matrix multiplication.

Termination Criteria

As mentioned earlier, we need to decide how to terminate an effort to align two $\beta$-strands. It is not aways possible to align two $\beta$-strands. The most common problem is that the two atoms we wish align are too far apart and no matter how we alter the backbone bonds between them they can never be parallel to each other. Less common, but still a problem is that they could be too close, and no assignment of torsion angles could cause the intervening turn region to bend sharply enough for the two atoms to interact.

To account for these situations where no solution can be found, we have two simple requirements as we iteratively attempt to set $${ \mathbf{T} ( \boldsymbol{\tau_{i+1}}) \doteqdot \mathnormal{G_\zeta}}$$. We require $${ \Vert\mathbf{T} ( \boldsymbol{\tau_{i+1}}) - \mathnormal{G_\zeta} \Vert < \Vert\mathbf{T} ( \boldsymbol{\tau_{i}}) - \mathnormal{G_\zeta} \Vert}$$, and we set some maximum number of iterations, usually ten. If either of these requirements is violated, we do not tweak the protein and continue with our normal search techniques.

Implementation Details for Tweaking

There are a number of implementation details that one must consider before the above techniques can actively be used to align $\beta$-strands into $\beta$-sheets. Here, we will highlight some of the issues that we need to deal with.

A critical step in each iteration of solving our program (5.20) is the construction of the $\mathbf{QR}$-factorization of our derivative matrix $\mathbf{A}$. We have chosen to use a freely available matrix solver called Meschach [72]. Using Meschach provides numerous advantages. First, the program is in ANSI C, so it is easily folded into our existing code. It was written by experts in numerical analysis, so the issues and considerations inherent in finding numerical solutions using computers have been considered. And lastly, one of the authors is a colleague at the University of Iowa and is accessible to answer questions and provide support.

Using Meschach makes our solution of these problems fairly simple. First, we load our matrices and vectors into the appropriate data structures. Next, we find the $\mathbf{QR}$-factorization using library functions, utilize the upper and lower triangular solver library functions, and finally use matrix multiply to return our solution to the iteration. All memory is dynamically allocated and de-allocated, so we can solve matrices of varying size without having to worry about a maximum size.

Tweaking is not attempted at every node of the search tree. When we enter a new node of the search tree, the torsion angles for the new amino acid are noted and we enter the subroutine attemptTweak (see Figure 5.6). This subroutine checks for potential tweaking events and alignments before actually calling strandsAlign which implements the algorithm and techniques discussed in this chapter.

AttemptTweak first determines if the current amino acid is the possible starting point for a $\beta$-strand. To do this, it will check the torsion angles of the newly added target amino acid to see if they fall within a certain range of standard $\beta$-strand torsion angles. If not, we will return to the main search tree and use the discrete torsion angles. Assuming the torsion angles are amenable, we then run through the previously folded protein portion to analyze the current secondary structure. When doing this, we are recognizing $\beta$-strands, $\alpha$-helices, and turn regions.

We then work backward from our target amino acid looking for a partner amino acid that is part of a $\beta$-strand. We require at least three intervening turn amino acids before this partner, so that enough flexibility exists in order to align the target and its partner. For our implementation, an amino acid is tweakable so long as it is not currently involved in a $\beta$-sheet or an $\alpha$-helix, and it has not been tweaked in a previously successful effort to build a $\beta$-sheet. This means that tweaking will not destroy existing secondary structure, but will work around it to create additional secondary structure.

As we search backward to the beginning of the protein for potential partner amino acids, once a potential partner is found, we attempt to create both a parallel and anti-parallel $\beta$-sheet. If an alignment is found, we increment our depth and move onto the next level in our search tree. If none is found, or we return after searching a branch of the tree, we continue to search for partners until we reach the first amino acid. Then we leave tweaking and continue on with the regular search algorithm.

As we mentioned earlier, we search the search space in a depth-first manner using a parameter called optimism (see Section 3.4). We also use optimism in our tweaking algorithm. When optimism is low we restrict how often we will attempt to perform a tweak. For instance, if amino acids are far apart or if the convergence is slow, we may terminate the tweak attempt before being certain if a tweak is possible. This meshes with our efforts in restricting the scope of the search until we have a higher level of optimism.

Figure 5.6: Algorithm: procedure $attemptTweak(protein,i)$.

Tweaking and Nagging

Tweaking must also work within our nagging infrastructure. If a child processor nags a parent and the parent's current state includes a successful tweak, the child must receive that information and rebuild the tweak as well. In the nagging case, the child will never backtrack over the tweak sent it by its parent, so it just sets all the angles the parent sends it, notes the angles which are the result of tweaks, and continues searching deeper in the search tree. This is possible because tweaking is a deterministic process. Both the child and parent processor will build the $\beta$-sheet in the same way.

Similarly, we have implemented a feature that allows us to rewind the search to a previous state if the root processor were to fail or is interrupted in its search. This is accomplished via a checkpoint file. The checkpoint file contains an output of all solutions found and an occasional output of the search state after some number of nodes. This output can then be read back in and the search state is updated to the latest search state in the checkpoint file. Tweak information must be stored and rebuilt as well. Unlike in the message passing for nagging, we have to rewind the search state when performing this type of restart. Therefore, we take great pains to properly travel through the subroutines setting all necessary counters and markers. This allows the search to properly backtrack and restart just as if it had never been halted in the first place.

Results from tweaking and examples of tweaked proteins will be given in Chapter 7.