Interestingly, although the method formally allows for fuzzy memberships, the optimal assignment of atoms to domains is always unique in practice, thus obtaining an exact partitioning of atoms into domains. For example, consider a hypothetical 3-atom system with the distance-deviation matrix(14)which has the optimal solution for M = 2 subunits:(15)Here atom 1 is placed in the first and atoms 2 and 3 are placed in the second domain. As described in Section “Cluster membership probability” the elements of the membership matrix either converge to one or to zero, i.e. the atoms tilt over to the domain that produces the smallest clustering error when including this atom. In other words, the clustering error is minimal when each atom is fully assigned to the subunit it belongs to most.

Now consider the case,(16)that has the solution for (17)with total error . In contrast the fuzzy solution(18)has a higher total error of . Finally, consider the pathological case of an off-diagonally uniform distance matrix which represents entirely uncorrelated motion(19)which may be found for gas particles. In this case the solution for is degenerated:(20)In this case the total error of was found for all . According to the uniformity of this degenerate case is never found in practice for macromolecules. Even nearly unstructured proteins will have some structure in because of their bonding topology and minor deviations in from the uniform case will cause the atoms to be uniquely assigned to one domain such that the error is minimum.

We conclude that no fuzzy memberships are found in macromolecules. Note, however, that the introduction of a fuzzy membership was still essential, because using this formulation we could express the optimization problem as a continuous quadratic optimization problem. The solution to this kind of problem is much easier than the solution to the integer optimization problem emerging by the priori assumption that the memberships must be integer values.

As a first example the optimization method was applied to Ala in order to demonstrate that the method can identify meaningful domains. Some of the resulting coarse-grain structures and the clustering error for are shown in Figure 1. The sub-structures are approximately equally sized and represent the optimal partitioning for a given number of domains. As the number of domains is increased the size of the domains diminishes. For , the method successfully identifies the domains that are nearly rigid due to bond angle, angle, improper dihedral and -angle constraints: There are 4 domains containing the 4 peptide planes including the first but excluding the second plus the side chain. The remaining two domains contain the N-terminal and the C-terminal (see Figure 1). The small remaining clustering error reflects the vibrations still allowed within the domains, mainly due to flexibility in the improper and -dihedrals. For the method clusters the system into domains containing one backbone atom each along with the one side chain atom connected to it. Finally the method is shown to be consistent in the limit, because for every atom is placed into a single domain (), and the error is zero (not shown as structure).

In order to study a more complex system, the method was applied to the MR121-GSGSW peptide. Figure 2 shows a series of molecular coarse-grain structures for selected numbers of domains . This series shows clearly how the flexible parts of the molecule subdivide into finer domains. separates the GSGS chain and the chromophores, splits also the GS domains. Using more domains accounts for smaller decrease of the error until the system is split into individual residues.

The corresponding distance-deviation matrix is shown in Figure 3 (top). It is structured into blocks along the diagonal (values are close to zero), that represent the almost rigid regions of the molecule. The values on the diagonal are zero () (blue), while values far from the diagonal are large (red). To identify rigid domains within the peptide we have employed the quadratic optimization method for . The convergence of the method depends on the size of the molecular system, the number of domains chosen and the initial conditions. For six domains in the artificial peptide (81 atoms) it converged within iterations, and a few on a standard desktop computer.

The resulting membership matrix, , and the corresponding coarse-grain structure are shown in Figure 3 (bottom). The colors show the assignment of atoms to domains. The elements of the membership matrix either converge to one or to zero, i.e. the atoms tilt over to the domain that produces the smallest clustering error when including this atom.

In order to test the optimality of the results, we have repeatedly solved the clustering problem for the MR121-GSGSW peptide using different initial conditions: (i) for given , each atom is assigned a random membership to each domain, and then normalized so that ; (ii) only is using a random initial condition while the solutions for are found by successive restart from the previous solution with the atoms of the largest-error domain split into the two new domains by an initial assignment of . Figure 4 shows a comparison for the clustering error for both cases, with ten realizations for the random initial condition plotted. The results are identical independent of the initial condition for small , which suggests that these solutions are likely to be globally optimal. For large the solution based on random initial condition gets trapped in different but only slightly suboptimal local minima, while the successive restart solution is monotonically improving for increasing . In all cases studied, the heuristic successive restart scheme possesses a useful monotonicity property, and performs better than optimization of random guesses.

The transport protein transthyretin (TTR) is primarily synthesized in liver, choroid plexus, and the retina. The primary function is the transport of thyroxine and retinol binding protein (RBP). Both molecules can bind to the homo-tetrameric structure of TTR, which is found at a physiological pH of . In contrast the kDa dimer structure is observed at pH and titration of sodium dextransulfate (SDS). It has two identical 127-amino-acid monomers (A - blue) and (B - green) (see Figure 5) with an extensive -sheet structure that form -sandwiches [43]. The interactions between the two monomers involve electrostatic and hydrophobic forces.

Transthyretin is one of the human proteins known to be associated with local amyloidosis. Amyloid fibrils are the polymerized form of the protein, their internal structure mainly consists of cross -sheets, arranged perpendicular to the long axis of the fibrils [44]. Both point variants of TTR and the native protein are known to deposit as amyloid fibrils in the extra-cellular region, where they cause neurodegeneration and organ failure (for reviews on amyloidosis see [45], [46]). Transthyretin is known to be associated with the amyloid diseases senile systemic amyloidosis (SSA), familial amyloid polyneuropathy (FAP), and familial amyloid cardiomyopathy (FAC). Other known amyloidogenic diseases are for e.g. Alzheimer's disease, type 2 diabetes and the transmissible spongiform encephalopathies which are characterized by proteinaceous deposits in the affected relevant organs.

Transthyretin aggregation to amyloid fibers has been the subject of many studies [43]–[48], however the molecular mechanisms are still not completely understood. Structural modifications and their effect on conformational stability were studied by structural and computational analyses [49] and experimentally by urea and temperature induced unfolding [50], [51]. It is proposed that amyloidogenicity of TTR is associated with anomalous structures that favour oligomer and fibril formation. The structures are assumed to be the product of complex dissociation via destabilisation [52] and subsequent unfolding and folding of the protein [53]. It could be verified that prior fibril formation the homo-tetramer dissociates into two dimers [54]. Whether the dimers need to dissociate into monomers before fibrillation can occur is still unclear. However it is assumed that dimer dissociation is the result of a mechanism called “edge exposure”, where the displacement of residues (inner -strand) and residues (outer -strand) flattens the dimer structure [55], [56].

To date, a large number of TTR variants have been associated with amyloid formation [57]. Here we study the structural rigidity of the wild-type and two variants commonly found, where the leucine of residue 58 in the dimers is replaced by arginine or histidine (TTR-58Arg and TTR-58His) and investigate their possible role in destabilisation and dissociation of the dimer structure. Both variants are known to be amyloidogenic, however the phenotypic difference of FAP between the 58His and 58Arg mutations suggest differences in the secretion efficiency or aggregation characteristics of the TTR variants [58].

In Figure 5 we show the clustering error for increasing number of domains and the coarse-grain structures for domains. As expected for the atoms from each monomer are placed in separate but symmetric domains. This separation is maintained for larger values of . For the domains found in the two monomers are nearly, but not perfectly symmetrical, as a result of limited statistical accuracy of the molecular dynamics trajectory. For the algorithm identifies two -sheets (b) and (d) in the dimer and two structures (a) and (c) including -sheets and the -helices. At the structures (a) and (c) are split into one block containing two -strands (f) and (g) and one block containing two -strands and the -helix (e) and (h). For the two loops in the outer region containing two short -strands are found to be separate domains (i) and (j).

To demonstrate the applicability to experimental data we used the method to partition molecular structures of transthyretin obtained by x-ray crystallography. Besides the wild type structure (PDB code 1DVQ), which was also used in the molecular dynamics simulation, five related structures of transthyretin complexed with resveratrol, diclofenac, flurbiprofen, DDBF, oFLU, and PHENOX (PDB codes 1DVS, 1DVT, 1DVU, 1DVX, 1DVY, 1DVZ) [39] have been used to generate the distance deviation matrix. Due to in sequence mutations in the structures 1DVX and 1DVZ we cleaned up the structure files to leave only comparable -carbon atoms in all six pdb-files. In the 1DVX file we removed residue 9+127 BLEU, 110+228 BSER and 113+231 BTHR, while in 1DVZ 7+124 BLYS and 10+127 BLEU where removed. We note that these six crystallographic x-ray structures correspond to different chemical or crystallographic states, so that the structural differences between them are not expected to be idendical to the structural differences within the Boltzmann-weighted ensemble of a solvated TTR in a single chemical state. Nevertheless, it is expected that the differences in the crystallographic realizations are sensitive to the molecule's instrinsic flexibility, so that a comparison between the simulation-based and X-ray-based results is interesting. The clustering result obtained from the x-ray structures and clustering error with coarse-grain structures for are shown in Figure 6. The clustering of x-ray structures yield similar coarse-graining as obtained by the clustering of molecular dynamics data. The dimer is separated into two monomers (A) and (B) for , while for the -sandwich (b)+(c)+(g)+(f) and -helical (e)+(h) structures within the dimers are identified (compare top right structure generated from MD data). However, for increasing number of domains () the clusterings are different. Note that the clustering error found for the X-ray structures is much smaller, indicating that the crystallographic realizations are much more similar to each other than the structures accessible to the dynamics of TTR in a solvent simulation, likely owing to crystal lattice constraints.

The distance-deviation matrices in Figure 7 show two large blocks along the diagonal (blue) that indicate that the internal rigidity within each of the two associated monomers is much larger than the rigidity between the monomers. The off-diagonal regions in the matrices (yellow-green-red) represent the inter-monomeric rigidity and are related to the stability or binding strength between the monomers. Large values in the matrix (red) indicate low stability, while small values (blue) are related to high stability of the dimer.

The structural modification induced by the amyloidogenic variants (TTR-58Arg and TTR-58His) contribute to an de/increase in rigidity in some regions of the structure (see increasing red regions in the variants compared to the native TTR in Figure 7), which lead to local de/stabilisation of the dimer. The overall stability of the dimer is directly related to the difference . Here, the binding strength for TTR-58Arg is increased and for TTR-58His is decreased compared to the binding strength of the wild-type . The decreased stability for TTR-58His variant compared to the wild type protein is supported by urea and thermal induced unfolding experiments [51] and computational studies that are based on an energy functions derived from non-redundant x-ray structures [49].

However in addition to the overall stability, increased atomic motion of specific regions in the dimer may influence the stability of the dimer and favor transient dissociation. The local flexibility/rigidity of atoms is reflected by , i.e. the mean row value of . The method is thus able to determine the relevant substructures that may cause destabilisation by taking the row average of the distance deviation matrix (see Figure 7). The peaks indicate residues that have increased distance deviation with respect to all other residues, i.e. the most flexible regions in the dimer. In Figure 7 (right) the structures are color coded according to the row average of , the mutated residue 58 is colored purple. In agreement with [51] the results indicate that compared to the wild-type protein the 58His and 58Arg variants are mainly destabilized at the monomer-monomer interface. In comparison to the wild-type TTR (Figure 7 top), it is clearly seen that the 58His variant increases the total distance deviation between the monomers (see average value of the row mean) and the peak values at residues . In contrast the 58Arg variant of transthyretin decreases the mean distance deviation, while the peak values at residues are still increased compared to the wild-type TTR. Because both variants are known to be amylogenic [58], we conclude that destabilisation is not only determined by the overall stability, but also by specific regions that cause local destabilisation of the protein that may lead to transient dissociation into monomers. The method is thus able to provide information about regions of TTR that are destabilized in disease causing variants.

The existence of semi-rigid domains and their relative dynamics are essential for the functionality of large macromolecular machines. Here, we analyse the dynamics of the GroEL-GroES chaperone complex (see Figure 8), which contains 8,015 residues (72,716 atoms). The complex ensures the proper folding of many proteins [59] and avoids non-native protein aggregation. GroEL is a tetradecameric protein of 14 identical domains arranged in a cis and trans heptameric-ring. GroES is dome shaped in either un-/bound configuration and contains seven identical domains assembled as a heptamer ring.

Computational studies have provided important insights into the allosteric mechanism of the chaperonin GroEL-GroES. Protein folding within the complex involves binding, encapsulation, and release of the substrate protein [60], [61]. During the GroEL-GroES cycle the GroEL binds a mis-/unfolded protein at its apical (A) domain (see Figure 9). The binding is caused by electrostatic and hydrophobic interactions between the exposed hydrophobic residues of the substrate protein and those of the apical domain. The equatorial domain (E) plays the major role in the overall chaperonin activity. It binds and hydrolyzes ATP. The intermediate domain (I) serves as a functional bridge between the apical and equatorial domains. After ATP binding to every cis-subunit, GroEL is bound to the cofactor GroES. During the GroES binding large conformational changes at the apical domain of the cis-ring cause upwards and outwards movement of the apical GroEL domains, thereby increasing the size of the central cavity and forming a dome-shaped chamber [59], [61]. By this conformational change the substrate protein is captured inside the cavity, where it will be able to undergo conformational changes toward the folded state. During ATPs hydrolysis in the cis-ring, ATP molecules are transferred to the trans-ring, which drives the release of the GroES cap and the substrate protein.

Due to the size of the system only a short molecular dynamics trajectory with duration of was produced, which is certainly not converged, but can nevertheless be used for a performance test (see Section on “Computational Performance”). In Figure 8 we show the clustering error for and the most informative structures. We note that due to insufficient statistical information in the MD simulations (under-sampling) the optimization results in disconnected domains (fragmentation) for small (see Figure 8). Nevertheless, the clustering with reveals the ring structure of the GroEL into two halves and finds the GroES as a third domain revealing the essential elements necessary to represent the conformational change caused by the complex formation of GroEL with GroES. For a large number of domains, e.g. for , the method clearly detects functional domains that can be directly related to the heptamer-ring structure of the chaperone. All shown results (, ) are equally “correct”, but reveal different levels of detail. The fragmentation of domains for small may be reduced by using longer molecular dynamics trajectories. Thus, the method is applicable to large molecular complexes with modest requirements of computation time, but as it is data-based the results are sensitive to the quality of the data.

To study GroEL in more detail, we have further performed domain identification on the subunit of the cis and trans heptameric-rings of GroEL. The distance deviation matrix was generated by averaging over the data of seven identical subunits for the cis and trans ring respectively. This averaging enhances the statistics of the molecular dynamics data. The domains found in the optimization (see Figure 9) are in good agreement with the functional domains in the GroEL subunit [62]. For the major three domains (apical, intermediate and equatorial) are found in either cis- or trans-ring subunits. The distance deviation matrix for and the identified domain boundaries are shown in Figure 10. The average row value and the color coding for the three domains is shown on the right. These coarse-grain structures are in good agreement to those used in rigid clusters models [63] or Markov models [64]. For larger number of domains, for e.g. , the method identifies two domains in the apical, intermediate and equatorial region respectively. For three domains in the apical and intermediate region respectively and four domains in the equatorial region were found.