Statistical Shape Methodology for the Analysis of Helices

Consider a helix in three-dimensional space along which a sequence of equally spaced points is observed, subject to statistical noise. For data coming from a single helix, a two-stage algorithm based on a profile likelihood is developed to compute the maximum likelihood estimate of the helix parameters. Statistical properties of the estimator are studied and comparisons are made to other estimators found in the literature. Next a likelihood ratio test is developed to test if there is a change point in the helix, splitting the data into two sub-helices. The shapes of protein α-helices are used to illustrate the methodology.

• Curved: There is a slow but steady change of the direction of the helix.
This can happen over a large part or even all of the helix.
• Straight: There is no change in the overall direction of the helix. Mardia et al. (1999) was another early paper on the estimation of curvature and torsion for a helix.
The study of helix structure falls within the scope of statistical shape analysis. Shape analysis deals with geometric objects in Euclidean space and is concerned with the properties that remain invariant under a group of transformations. There is now a rich collection of statistical tools for shape data (e.g. Dryden and Mardia, 2016). The simplest type of object consists of a collection of points or landmarks, and the most important groups are the similarity transformations (location, scale and rotation) and the rigid body transformations (location and rotation). For this paper an object consists of a set of points in R 3 that lie on or near a helix and the relevant group is the rigid body transformations.
A helix in three dimensions is defined by the function where r and 2πc denote the radius and pitch. In this representation the helix is traversed at constant speed as a function of an independent variable t, which can be regarded either as arc length after projecting the helix onto the plane of the first two coordinates, or as a sort of time, even though the helix exists as a static object. In addition, the position of the helix in R 3 can be altered by a rigid body motion. In a statistical helix, regular measurements are available along the helix subject to statistical error. There are two main purposes for the present paper. First we revisit the estimation problem for a statistical helix and show that maximum likelihood estimation can be reduced to an optimized least squares problem under certain assumptions. Secondly, we present a new global approach to the change point problem and investigate the use of feature statistics to highlight the character of a change point. S10

M.F. Alfahad et al.
In more detail, the paper is organized as follows. Section 2 introduces the statistical helix model (similar in spirit to Mardia et al. (2018), but more geometrically explicit). Section 3 gives the details behind the optimized least squares algorithm. This algorithm needs an initial estimate of the helix axis to start the iterations, and there are a large number of methods in the literature such as Rotfit described by Christopher et al. (1996) and HELFIT by Enkhbayar et al. (2008). Two methods to estimate the initial axis are compared in Section 4: Rotfit and a new method based on modified principal components. Section 6 presents a new likelihood ratio test for the presence of a global change point, both where the time index of the potential change point is known and where it needs to be estimated. We give this procedure the name ChangePoint-Detector. A bootstrap procedure is proposed to assess statistical significance. This procedure is investigated on simulated data in Section 7. In Section 8 it is applied to several protein examples and compared to Kink-Detector (Mardia et al., 2018).

The Statistical Helix Model
If an arbitrary rotation and location are incorporated in Eq. 1.1, then the mathematical helix at a regularly spaced set of "times" takes the form f (t i ) = r cos(t i )u + r sin(t i )v + ct i w + b, i = n 1 , n 1 + 1, . . . , n 2 − 1, n 2 .
(2.1) Thus the number of points on the helix is n = n 2 − n 1 + 1. General indices n 1 < n 2 are allowed for the start and finish points in order to facilitate the parameterization when a change point is present; see Section 6. Here • Γ = u v w is a 3 × 3 orthogonal matrix whose three columns define the orientation of the helix. In particular the vector w defines the helix axis, and the vectors u and v define the plane normal to the helix axis.
• r > 0 defines the helix radius, • 2πc > 0 defines the helix pitch, which is the vertical height of one turn of the helix, • b ∈ R 3 is an intercept, • t i = iβ is a sequence of regularly-spaced times at which the helix is observed, where β > 0 defines the turn angle of the helix in radians (i.e. the angle between two consecutive points on the helix).

S11
Helices can be regarded as right-handed or left-handed depending on whether det(Γ) = +1 or −1 respectively (e.g. Campbell and Farrell, 2014). For the purpose of this paper we largely restrict a helix to be right-handed.
A regular (statistical) helix is a collection of points or landmarks in three dimensions obtained from a mathematical helix by adding noise, where the error terms are assumed here to follow independent isotropic normal distributions. The adjective "regular" is used to distinguish model (2.2) from the change point helix model to be introduced in Section 4. It is also convenient to let so that g (t) − b denotes the projection of the centered true helix function f (t) − b onto the helix axis w . In particular, let denote the axis values at the data times t i . For the purposes of this paper, we shall treat the turn angle β as known. It is well-known in the protein structure literature that the turn angle β along the helix can be treated as having a constant value very close to β = 100 • ; e.g. Dickerson and Geis (1969) give a value of 1.75 radians = 100.3 • . Recent confirmation can be obtained from the detailed analysis of the maximum likelihood estimation of β carried out for 129 straight helices from crowdsourcing data in the web supplement (Web Figure 5(d)) to Mardia et al. (2018); for this data it was found that the mean estimate of β is 99.1 • with a standard error of 1.2 • .
All of the other parameters will be regarded as unknown and needing estimation. However, for the purposes of developing an estimation algorithm, we shall first treat the case where the axis w is known.
The parameters of a helix can be divided into two types. The registration parameters are the orthonormal vectors u, v and w , and the intercept vector b = (b 1 , b 2 , b 3 ) T . The shape parameters are the radius r and the pitch c. Consider a right-handed helix in Eq. 2.2 with orientation matrix Γ, so Γ T Γ = I and det(Γ) = +1. The following definitions impose further restrictions on Γ. The helix is said to be in is the identity matrix, so the three orientation vectors are given by the three standard coordinate directions in R 3 ; in particular w 0 = 0 0 1 T lies in the vertical direction. For a left-handed helix in canonical coordinates it is convenient to define Thus, looking at the horizontal plane from above in canonical or semicanonical coordinates, a right-handed helix winds around the plane in a counter-clockwise direction as t increases; a left-handed helix winds in a clockwise direction.
3 Parameter Estimation for a Regular Helix 3.1. Known vertical helix axis Start with the assumption that a righthanded data helix is in semi-canonical coordinates, so that w = w 0 = 0 0 1 T is known to be vertical. Then u and v take the form for some angle τ . In this case model (2.2) can be re-written as Statistical Shape Methodology...

S13
where α 1 = r cos τ, α 2 = r sin τ . The model in Eq. 3.2 can be viewed as a multivariate linear regression model with a three-dimensional response on n observations. The regression parameters are α 1 , α 2 , c, b. Since the error term is isotropic, the model can also be represented as a multiple regression model with 3n scalar responses, after stacking the 3 columns for the response. Further, maximum likelihood estimation reduces to least squares regression. The 3n × 6 design matrix is Write y i = (y i1 , y i2 , y i3 ) T . Then the least squares estimators take the form (3.4) in terms of the centered variables We can then deriveτ andr bŷ where atan2 is the two-argument arctan function, so that we have (α 1 ,α 2 ) = r(cosτ , sinτ ). Further, the estimated shift vectorb = [b 1 ,b 2 ,b 3 ] T is given byb where y 1 , y 2 and y 3 are the means of each coordinate. Finally, the residual sum of squares (RSS) is given by: whereŷ i is the i th fitted value of y i , a vector of dimension 3. The residual sum of squares depends on the choice of helix axis, denoted here by w 0 .
In the case of a left-handed helix, just change the sign of one of the columns in Eq. 3.1; e.g. let for some angle τ , with corresponding changes to the subsequent algebra. If it is unknown whether or not the helix is right-or left-handed, fit both models and choose the model with the smaller residual sum of squares. Unless σ 2 is extremely large, the correct choice will be obvious. In addition, if the estimate of the pitch parameter c in Eq. 3.4 is negative, then it is necessary to change the sign of the helix axis w (plus the sign of either u or v to preserve the sign of det(Γ)).

Known helix axis in general position
Next let the axis w of a righthanded helix be a known unit vector, but not necessarily vertical. Let G = G(w ) be a 3 × 3 rotation matrix whose third column equals w . Let z i = G T y i , denote the rotated data, so that the known helix axis for the {z i } is vertical. Then the estimation procedure of the previous section can be applied to the {z i }.
The plane spanned by the first two columns of G is determined by w , but not the two columns themselves. A rotation of this plane about the w axis corresponds to changing the meaning of the angle τ in Eq. 3.2.
After fitting the helix by least squares, the quality of the fit can be summarized by the residual sum of squares, denoted RSS(w ), say. The value of RSS(w ) does not depend on the indeterminacy in the meaning of τ .
3.3. Unknown helix axis If w is unknown, an iterative method based on profile likelihood can be used to find the maximum likelihood estimators. The procedure works as follows: (a) Start with an initial estimate w init , say. Two possibilities are suggested in the next section.

S15
(b) Given w , the maximized likelihood with respect to the remaining parameters (known as the "profile likelihood"), is a monotone decreasing function of RSS(w ). A nonlinear optimization algorithm, e.g. the routine nlm in R (R Core Team, 2014), can be used to numerically minimize RSS(w ) over w lying on the unit sphere in R 3 . There is no mathematical guarantee that RSS(w ) possesses a unique local minimum. Hence it is important to choose a good starting point. For all the examples in the paper, convergence has never been as issue.
For the optimization in (b) it is helpful to rotate the initial estimate w init to the north pole, w 0 = 0 0 1 T , and to represent w using an As (p 1 , p 2 ) ranges through R 2 , w covers the unit sphere minus the south pole, 0 0 −1 T . In practice the minimizing value of w will usually be very close to the initial estimate w 0 . The resulting MLEs for all the parameters can be termed the "Optimized least squares" (OptLS) estimates. An estimate of the error variance iŝ where in the denominator we subtract 8 degrees of freedom from 3n since a regular helix model contains 8 regression parameters (6 for the linear regression model, plus 2 for the helix axis).

Initial Estimate of the Helix Axis
Several methods have been established in the literature to estimate a helix axis; see for example, Christopher et al. (1996) and Wilman (2014a). Here we limit discussion to two methods for getting initial estimates: Rotfit S16 M.F. Alfahad et al.
(described in Christopher et al., 1996) and a new method based on modified least squares.
The principle behind Rotfit is easy to describe. Starting from an n × 3 mathematical helix Y , let Y −1 denote Y without its first row, and Y −n denote Y without its last row. Then Y −1 can be mapped onto Y −n by a shift and rotation. Further the fixed axis of the rotation matrix is the desired axis w .
When working with a statistical helix, the rotation matrix can be fitted using Procrustes analysis (e.g., Mardia et al., 1979 p. 416) where M and N are 3 × 3 orthogonal matrices and L is a diagonal matrix with positive entries. Then the Procrustes rotation matrix can be estimated by R = MN T . Further, the fixed axis of R is the eigenvector with eigenvalue 1, and can be found by a spectral decomposition (the other two eigenvalues are complex).
The modified least squares method can be described as follows. Starting from Y , construct the increments and combine them into an (n − 2) × 3 matrix D.
In a mathematical helix, E will have a zero eigenvalue with eigenvector given by w . In the statistical case, there will be one small eigenvalue, and two larger approximately equal eigenvalues.
Recall that if w is an eigenvector, so is −w . Hence we need to specify its sign. That is, we need to choose the sign of w so that the fitted value of pitch parameter c in the helix model (2.1)-(2.2) will be positive. This task is straightforward when the level of noise σ 2 is not excessive. Just ensure the sign of w is chosen so that the difference between the endpoints, after projection onto the helix axis, is positive, w T y n − w T y 1 > 0.
In practice, it does not matter which of these two initial estimates is used. However, simulations from the model (2.2) indicate that Rotfit is generally more accurate.

Assumptions About the Turn Angle β
The turn angle β is a key parameter in the statistical helix model, and there are several ways in which it can be treated.
Model A: the helix turn angle β = 100 • is known exactly. This is the choice made in this paper; see Section 2.

S17
However, as emphasized by a referee, it is also of interest to consider what happens when this assumption is violated. There are two natural possibilities: Model B: the turn angle β is constant within a helix, but is not known exactly. One course of action is to include it as one of the parameters to be estimated in the likelihood for a single helix, as in Mardia et al. (2018). Another is to carry out a sensitivity analysis to assess the effect on the conclusions of varying β. We carried out a small sensitivity analysis on simulated data and found that varying β in the range 98 o − 102 • had a negligible effect on the statistical analysis. This range has been chosen to be roughly 100 o ± 2s where s = 1.2 • is the standard error reported in Section 2 that was found by Mardia et al. (2018).

Model C:
A more severe violation is to allow the turn angle to vary along the helix; that is, the incremental turn angle between x j and x j+1 , denoted now by β j , varies with j, j = 1, . . . , n − 1, where n is the number of landmarks in the helix. The simplest explicit model is to assume that the β j are i.i.d. samples from N (β 0 , τ 2 β ), say, with β 0 = 100 • .
For illustration suppose τ β = 4 • . This choice is motivated from Web Figure 5(d) in the web supplement to Mardia et al. (2018), where the 129 straight helices have lengths ranging from 16 atoms to about 40 atoms with most helices having 20-25 atoms. Since the standard deviation of the 129 values for β is about 1 • , and since each of these values stems from a data set with n ≥ 16, one can conclude that the standard deviation of the distribution of the β j is at least τ β = 4 • .
Next assume we have a helix consisting of n = 2m + 1 atoms, where the central atom is perfectly aligned with the model. The outermost atoms have cumulative turn angles For example, this means that the outermost atoms in a helix consisting of 19 = 2×9+1 atoms will have turn angles which are distributed with a standard deviation of √ 9 × 4 o = 12 • around the position assumed by a model with a constant turn angle. Hence errors of 10 -15% in the cumulative turn angle can easily occur under Model C. S18 M.F. Alfahad et al.
Model C was explored in the web supplement to Mardia et al. (2018), and was found to be unsupported by the data in the protein application. The view in biochemistry is that a constant incremental turn angle is more appropriate than Model C. In addition the external knowledge about β has a small enough standard error to justify using Model A instead of Model B.

The Change Point Model
We have already used the term regular helix model to describe the statistical helix in Eq. 2.2 in which all the data are modelled by a single helix. In this section we introduce and investigate a change point helix model in which there is a change point along the helix. That is, for some value of k, n 1 < k < n 2 , we suppose that the points y i , n 1 ≤ i ≤ k lie on one statistical helix, and the remaining points y i , k + 1 ≤ i ≤ n 2 lie on another helix. All the regression parameters of the two helices are allowed to be different from one another, but they are assumed to have a common value of the error variance σ 2 . Note that the change point model does not require the continuity of the curve consisting of the two helix pieces at the change point t = k + 1 2 . Let m denote the smallest number of points along a helix we are willing to countenance for estimation purposes. Then the admissible values of k are Since a regular helix model contains 8 regression parameters, the smallest feasible choice of m is m = 3. However, it is reasonable to limit attention to larger values of m in order to avoid very short helices. Following Mardia et al. (2018), we take m = 6 in this paper.
6.1. Known change point index k Suppose the index k of a possible change point is known, and consider testing for the presence of a change point. The null hypothesis is H 0 : the data follow a regular helix model, and the alternative hypothesis is H 1 : the data follow a change point helix model.
Just as in the analysis of variance, the likelihood ratio test can be recast in terms of an F statistic. Let RSS (0) denote the optimized residual sum of squares under the null model after fitting a single helix to the whole data set by OptLS with d (0) = 3n − 8 degrees of freedom. Similarly let RSS ( ) (k) denote the optimized residual sum of squares after fitting a helix to the data points i = n 1 , . . . , k for = 1 and i = k + 1, . . . , n 2 for = 2, respectively. Then the likelihood ratio test statistic can be recast as the F -statistic (6.1) Statistical Shape Methodology...

S19
Here SS is the reduction in the residual sum of squares after fitting the alternative model with degrees of freedom d ( If the error variance σ 2 is small, we expect that the dependence on the parameters in the helix model (2.2), and its change point counterpart, can be linearized through a Taylor series expansion, so that the model can be approximated by a standard linear model with normal errors. Hence from the standard ANOVA decomposition, we expect that when the change point position is known, the F k -statistic will approximately follow the F (d (B) , d (W ) ) distribution, as suggested by the notation in Eq. 6.1. To test this expectation, a simulation study was carried out with n = 30, r = 2.3, c= 5.4 2π and σ 2 = 0.05. Then 10, 000 helices were simulated from the null distribution and the statistic F k was computed to look for a change point at k = 15. A Q-Q plot comparing the simulated F -statistics to the F (8, 74) distribution is given in Fig. 1. Note that there is close agreement between the two distributions except in the upper tail, where the simulated test statistic has a shorter tail than the F -distribution.
6.2. Unknown change point index k In practice the location of the change point is generally unknown. The likelihood ratio statistic is now given by maximizing F k over k, where m = 6. Since k is no longer fixed, the distribution of F max will be larger than the F (d (B) , d (W ) ) distribution. Hence standard statistical tables S20 M.F. Alfahad et al. Figure 1: The Q-Q plot of 10,000 simulated F k -statistics from Eq. 6.1 versus quantiles from the F-distribution for the example in Section 6.1 cannot be used for testing purposes. Instead the parametric bootstrap is used to assess significance. The parametric bootstrap works as follows. Under the null hypothesis, estimate the regression parameters and σ 2 . Using these values, simulate a large number n boot of new helices with normally distributed error terms. For each simulation evaluate the value of F * max , where * indicates a simulated bootstrap sample. Then the upper α tail of the simulated distribution, F * (α) max , say, gives the threshold of a statistical test of size α, e.g. with α = 0.05.
If the null hypothesis is rejected, then we can try to pinpoint the ways in which the two sub-helices differ from one another. Letk denote the value of k maximizing (6.3). The following notation is useful. LetΓ = ûvŵ denote the fitted orientation matrix under the null hypothesis, and letĝ ( ) (t) andf ( ) (t) denote the fitted axis and helix functions for each sub-helix.
ThenΓ Tĝ ( ) (t) andΓ Tf ( ) (t) denote the fitted axis and helix functions after rotation to canonical coordinates under the null hypothesis. Let denote the fitted axis position and fitted helix position, respectively, of a notional landmark at indexk + 1 2 under each of the sub-helix models, after rotation to canonical coordinates under the null hypothesis.
The following six features can be constructed to compare fitted subhelices. Note that they are not orthogonal to one another.

(Difference in helix axes.) Let
A 1 = 1 − cosθ whereŵ (1)Tŵ (2) = cosθ and where 0 ≤θ ≤ π is the angle between the two sub-helix axes. The statistic A 1 contains the same information asθ and measures the difference in helix axis for the two sub-helices, = 1, 2. If the two helices point in the same direction thenθ = 0 and cosθ = 1 so that A 1 = 0. Further, A 1 increases as the directions get further apart. The angleθ is the most important feature in practice when there is a change point.

(Shift parameter.) Let
3 ) 2 denote the squared difference in the fitted axis positions at indexk + 1 2 , as measured along the fitted helix axisŵ under the null model in canonical coordinates.

(Offset parameter.) Let
2 ) 2 denote the squared Euclidean distance between the fitted axis positions at indexk + 1 2 , after projection onto the (û,v )-plane under the null model, which is a horizontal plane in canonical coordinates.

(Spin parameter). Let
1 ), denotes the angle, treated as a number in [−π, π), between the two fitted helix angles, after projection onto the (û −v )-plane under the null model. This feature measures spin at the change point between the two sub-helices.

(Change in pitch.) Let
Thus we have six features: difference in helix axes, shift parameter, offset parameter, spin parameter, change in radius, and change in pitch. In each case critical values for these features can be simulated using the parametric bootstrap. For all quantities (A 1 , A 2 , A 3 , A 4 , A 5 , A 6 ), an upper critical value is used. In each bootstrap sample, the values of A 1 , . . . , A 6 are computed using the bootstrap value of k = k * maximizing the bootstrap version of the F -statistic given by Eq. 6.3.
A simple simulation study was carried out to see how the distribution of F max depends on σ 2 and n using parameters typical for protein α-helices. On fixing n = 30 and varying σ 2 , with n boot = 1000, we found that the bootstrap threshold values for F * max , in Table 1, are very similar for all σ 2 , which implies the distribution of the F statistic does not depend heavily on σ 2 . On the other hand, if we fix σ 2 = 0.05 and alter n, then the threshold value does change somewhat; that is, the distribution of the F max statistic does depend on the number of landmarks, just as the usual F distribution does.
For each data helix we can compare the computed F statistic with the threshold, F * (α) max . If F < F * (α) max , we conclude that there is no evidence our helix has a change point; otherwise, if F ≥ F * (α) max , the point with the largest F statistic corresponds to the estimated index of the change point. We give our procedure for estimation and testing the name ChangePoint-Detector. Statistical Shape Methodology...

Simulation Study for ChangePoint-Detector
In this section we illustrate the behavior of ChangePoint-Detector on two simulated helices, a regular helix and a change point helix.
• Example 1. A regular helix was simulated with n = 27 landmarks and parameters that mimic a protein α-helix (i.e. r = 2.3, 2πc = 5.4, β = 2π 3.6 and σ 2 = 0.05). • Example 2. A change point helix with change point k = 12, was constructed by simulating a regular helix in the same way as Example 1, but then introducing a rotation of the helix axis by an angle θ = 0.3 radians = 17.1 • , centered at time index k + 1 2 = 12.5. Table 2 presents the results for the key parameters, the overall test statistic and the six features of interest. The p-values in parentheses are constructed using bootstrap sampling with n boot = 1000.
First consider Example 1, the regular helix. As expected, there is no reason to reject the null hypothesis. The test statistic F max is not significant, the estimates of σ 2 under the null hypothesisσ 2 = 0.039 and under the alternative hypothesis σ 2 p = 0.039, from Eq. 6.2, are approximately equal, and none of the features A 1 , . . . , A 6 are significant.
Next consider Example 2, the change point helix, for which the test statistic F max is highly significant. For this example, ChangePoint-Detector   Table 3 shows the bootstrap distribution of k * for n boot = 1000 simulations to help judge the accuracy ofk. In particular, k * equals the estimatek = 12 most of the time, but occasionally overshoots by 1 or 2. The estimated angle between the two sub-helicesθ = 0.32 radians is close to the true value. Further the p-value for A 1 confirms that the change in angle is significantly different from 0, but none of the other features (A 2 , A 3 , A 4 , A 5 , A 6 ) is significant. Figure 3 shows the fitted helices. Panel (a) shows a fitted single helix for Example 1. Panel (b) shows the two fitted sub-helices for Example 2 after estimating the change point. The change in helix axis of sizeθ = 18.4 • = 0.32 radians at k = 12 is visible.   Figure 3: The two sub-helices found by ChangePoint-Detector for the simulated regular helix (Example 1) and the simulated change point helix (Example 2). In each case the two sub-helices are plotted using a thin line and a thicker line, respectively. Also plotted are the two fitted axes, as a dotted line with small dots and a dotted line with larger dots, respectively 8 Data Examples In this section we look at nine α-helix protein structures from Mardia et al. (2018). Each helix comprises a sequence of C α atoms; there are 3.6 equally-spaced C α atoms per turn of the helix, i.e. one C α atom every β = 2π 3.6 radians (100 degrees). For the i th landmark (C α atom) we associate the time index t i = iβ, i = n 1 , . . . , n 2 , where all the helices start at n 1 = 1 and and for each helix n 2 = n denotes the total number of the landmarks.
Three other parameters in the helix model in the protein setting have ideal values: 2πc = 5.4Å (the vertical distance of one complete turn), radius r = 2.3Å, and the variance σ 2 = 0.056Å 2 (e.g., Mardia et al., 2018). One Angström (Å) equals 10 −10 m. For the analysis in this paper, we treat β as known, but include r, c and σ 2 as unknown parameters in our estimation and testing procedures, so that the data can "speak for themselves" as far as possible.
In the study of protein helices, it is of particular interest to identify and locate kinks. Typically the number of kinks in protein helices is either 0 or 1 (Mardia et al., 2018). Mardia et al. (2018) developed an algorithm called Kink-Detector to identify the presence of a kink and to estimate its location. The Kink-Detector method is based on a local moving window (6 landmarks each side of the possible kink).

M.F. Alfahad et al.
Following Mardia et al. (2018), we used their nine test helices here, labelled Helix 1 -Helix 9. Mardia et al. (2018) chose these helices partly because identification of any kinks was particularly challenging. Note that the conclusions from the Kink-Detector methodology were confirmed in a crowdsourcing assessment (Wilman et al., 2014b) by experts in the protein community.
Here we want to see if the kinks found by Mardia et al. (2018) coincide with any change points found by ChangePoint-Detector, or if the presence of a kink is unassociated with larger-scale structure in the helix.
The ChangePoint-Detector methodology was applied to each of the nine helices. The regular helix model H 0 consisting of a single helix and the change point helix model H 1 with one change point were fitted to the data. The F max statistic of Section 6.1 was computed to test H 0 vs. H 1 . Further the statistics A 1 to A 6 were computed to compare the two sub-helices under H 1 . The results are summarized in Tables 4-5. These tables show for each helix: (a) the error variance estimatesσ 2 andσ 2 p from Eq. 6.2 (b) the overall F max statistic; (c) the optimizing indexk for the change point; (d) the anglê θ in degrees between the two sub-helix axes; and (e) the features A 1 , . . . , A 6 ; and (f) the p-values for F max and A 1 , . . . , A 6 based on bootstrap sampling with n boot = 1000.  From Tables 4-5 we see that all nine helices have a highly significant value of F max (with p-value in each case less than 0.001), meaning that all helices are deemed to have a change point. Further, in general the change point seems to be due to a change in axis direction. For the most of the helices, the feature A 1 is significantly different from 0 with a p-value less than 0.001. Two mild exceptions are Helices 5 and 6 with p-values 0.008 and 0.049. Some further evidence in support of the change point helix model is given by the residual variances. In most casesσ 2 p is close to the theoretical value of σ 2 = 0.056. The main exceptions are Helix 2 withσ 2 p = 0.144 too large and Helix 8 withσ 2 p = 0.014 too small. Figure 4 presents Helix 8 and the fitted helices using ChangePoint-Detector. The axis has a change point witĥ θ = 9.6 • as the first sub-helix axis direction (bottom) is different than the second sub-helix axis direction (top). Figure 5 shows that the maximum F k occurs at k = 8.  Figure 4: The two sub-helices found by ChangePoint-Detector for Helix 8. The two sub-helices are plotted using a thin line and a thicker line, respectively. Also plotted are the two fitted axes, as a dotted line with small dots and a dotted line with larger dots