Dynamics and Thermodynamics of Transthyretin Association from Molecular Dynamics Simulations

Molecular dynamics simulations are used in this work to probe the structural stability and the dynamics of engineered mutants of transthyretin (TTR), i.e., the double mutant F87M/L110M (MT-TTR) and the triple mutant F87M/L110M/S117E (3M-TTR), in relation to wild-type. Free energy analysis from end-point simulations and statistical effective energy functions are used to analyze trajectories, revealing that mutations do not have major impact on protein structure but rather on protein association, shifting the equilibria towards dissociated species. The result is confirmed by the analysis of 3M-TTR which shows dissociation within the first 10 ns of the simulation, indicating that contacts are lost at the dimer-dimer interface, whereas dimers (formed by monomers which pair to form two extended β-sheets) appear fairly stable. Overall the simulations provide a detailed view of the dynamics and thermodynamics of wild-type and mutant transthyretins and a rationale of the observed effects.


TTR: Origin and Function.
Formerly known as prealbumin, transthyretin (TTR) is a 55 kDa globular oligomeric protein made up by four identical monomeric units (I-IV, see Figure 1(a)) each composed of 127 amino acid residues. Before being released into plasma, TTR is principally produced in liver, choroids plexus, and retina [1][2][3][4]. It is mainly a carrier protein that binds to retinol binding protein (RBP) to transport vitamin A in plasma and represents the leading transporter of thyroxine (T4) in cerebrospinal fluid (CSF) [5][6][7].
1.2. TTR: Structure. The 3D-conformation of all TTR-tetramers displays high degree of symmetry as shown at high resolution by Blake and coworkers [8] and Hörnberg et al. [9]. TTR is an overall -sheet protein with a smallhelix domain between strands E and F (see Figure 1(a)).
Its structure is highly ordered with a flexible N-terminal region that could not be resolved in many crystallographic studies. In each TTR monomer, about 45 % of residues are arranged in a sandwich immunoglobulin-like topology made up by two four antiparallel-strand -sheets; the inner sheet DAGH is opposed to the outer sheet CBEF. Its secondary structure also exhibits a short -strand portion A * (see Figure 1) which is antiparallel (folded back) to strand A through a -turn (i+5) and which is involved in dimer-dimer contact [10]. The dimer of monomers I and II (Figures 1(b) and 1(c)), also known as primary dimer, is the crystallographic asymmetric unit structurally symmetric to the dimer of monomers III-IV. It is stabilized by a network of six main chain hydrogen bond interactions established by three pairs of residues (A120(NH) (OC)Y114, T118(NH) (OC)Y116, Y116(NH) (OC)T118) at the HH' interface ( Figure 2). Monomer I inner sheet DAGH forms , and S117 (purple-blue)) are shown explicitly (stick representation) as well as the label of different -stranded regions schematised in Figure 2, PDB id 1F41 [9]. S46  A45  F44  P43  E42  W41   A29  V30  H31  V32  F33  R34  K35   L55  E54   G53   L12  M13  V14  L15  V16  L17  D18   A19   T123  V122  V121  A120  T119  T118  S117  Y116  S115   Y114   R104  Y105  T106  I107  A108  A109  L110  L111  S112   P113   A91  E92  V93  V94  F95  T96  A97   F   H90   E89  I73  E72  V71  K70  Y69   I68  G67   E   D74  F87   S115  Y116  S117  T118  T119  A120  V121  V122  T123   Y114   A97  T96  F95  V94  V93  E92  A91   H88   T75   H90  E89  H88  ). Letters A to H designate the -strands in one monomer while H' and F' belong to the other monomer of the same dimer. Red arrows indicate the connection between strands and show the loop regions and blue ones go from donor (NH) to acceptor (OC). Coloured residues (yellow) are the ones involved in the mutation and is the helical region between strands E and F. with its homologue H'G' A'D' in monomer II a kind of symmetric antiparallel pseudo-continuous eight-strandsheet centered on HH' (DAGH-H'G' A'D'), Figure 1(b) [10,11]. At the opposite side, the pseudo-continuity in outer sheets CBEF-F'E'B'C' (Figure 1(c)) is rather loose and only 4 backbone hydrogen bonds (F87(CO) (HN)T96 and E89(NH) (OC)V94) are established ( Figure 2). While dimers are mainly kept together by hydrogen bonds at monomer-monomer interface, tetramers become the dimers of dimers essentially stabilized by hydrophobic contacts established within neighbouring subunits at the dimer-dimer interface via AB-GH loop interactions involving mainly residues L17, A19, V20, L110, P113, T119, and V121.

TTR-Related
Amyloidosis. TTR-related amyloidosis can be inherited in the case of genetic mutations or can be nonhereditary when it is due to wild-type (WT-TTR) [12]. WT-TTR amyloidosis also termed senile systemic amyloidosis (SSA) is the prevalent form of TTR amyloidosis which principally affects the heart. It develops with ageing and requires TTR tetramer dissociation and partial unfolding. Mutation-induced TTR amyloidosis is instead associated with familial amyloid polyneuropathy (FAP) which impacts essentially on peripheral nervous system and familial amyloid cardiomyopathy (FAC) with broad implication on heart [7,12]. The interconnection between protein structure and stability and its capability to form amyloid fibrils have motivated several structural studies. To date, about 80 point mutations have been correlated to human inherited amyloidosis [13].

Designed Mutants.
WT-TTR displays a very stable tetrameric molecular assembly. It is established that its structure is resistant to dissociation at physiological concentrations within the pH range 5-7 [14]. Several amyloid fibril formation models have been proposed and reviewed [15][16][17]. However, there is almost a consensus that the common mechanism involves the dissociation of the native tetrameric TTR into unstable but folded monomers, followed by local unfolding of the latter into multiple nonnative amyloidogenic intermediate states that self-assemble in solution [2,6,[15][16][17][18]. Furthermore, evidence has been provided pointing to the dissociation of the native tetramer as the rate-limiting step towards aggregation [5,[19][20][21][22]. Therefore knowledge of the dynamics and thermodynamics of dissociation of TTR native structure is an important issue as tetramer dissociation, monomer misfolding, and self-assembly of amyloidogenic monomers into amyloid and other aggregate morphologies are known to be linked to several human degenerative diseases. Specifically, TTR represents one of the few examples whereby it was possible to devise a drug for TTR amyloidogenesis [23] that acts through the stabilization of a native structure, thus indicating that thermodynamics studies and modeling of tetramer dissociation are a very important target.
To mimic the pathological situation, previous sitedirected mutagenesis, namely, mutation F87M/L110M (MT-TTR), was carried out to promote the dissociation of TTR tetramer into monomers [5]. To further shift the equilibrium towards monomers to obtain a more homogeneous dissociated species, an additional amino acid replacement has been introduced in MT-TTR molecule (S117E) [24]. Indeed, the latter mutant containing three mutations (3M-TTR) proved to be markedly more prone to in vitro monomerization in comparison with the double mutant MT-TTR, whose tetrameric states could be recovered by Tafamidis binding, at variance with 3M-TTR [24].
The three mutations are not reported in the ClinVar database [25], although pathological mutations do occur in neighbouring residues.
1.5. TTR Amyloidosis Therapeutics. As mentioned previously, the dissociation of TTR tetramer is believed to be the initial step into its fibrillation pathway [26]. Therefore, many strategies to TTR amyloidosis prevention exploit its ability to bind small molecules in the T4 binding channel (see Figure 1(a)), mimicking its hormone binding capability, thereby producing kinetic stabilization of the tetramer [27,28]. Other clinical remedies to FAC and FAP amyloidosis often employ organs transplantation (e.g., liver), even though not all the affected organs can be transplanted (e.g., choroids plexus, where TTR is produced as well [29]). Furthermore, additional treatments against TTR amyloidosis have been experimented, among which there are gene therapy [30] and the resorption of amyloid deposits [31].
The aim of this work is to provide a rationale of the observed dynamics and thermodynamics of the wild-type and mutant transthyretin based on atomistic molecular dynamics simulations. The methods applied here show the relative stability of monomers, dimers, and tetramers and provide a description of both enthalpic and entropic contributions taking advantage of recently developed methods.

Free Energy from End-Point MD Simulation.
The free energy of a system, with respect to some reference state, can be in principle deduced from conformational samples obtained, e.g., in a molecular dynamics simulation.
For a well equilibrated system, the enthalpy can be estimated as an ensemble average of the potential energy over a set of conformations from the MD trajectory. The bottleneck in obtaining the Gibbs free energy is the proper modeling of the entropic part and in particular of the solvation entropic effects.
In order to obtain free energy estimates all molecular dynamics snapshots have been analyzed using the GBSA implicit solvent model as implemented in a home-written version of the software Bluues [32,33].
The Gibbs free energy in the context of implicit solvent MD simulation has been discussed in reviews [34][35][36][37]. We follow here our recent perspective on the issue [37]. The free energy of a system may be written as where is the solute potential, + ΔW is the potential of mean forces (with ΔW the solvation energy), B is the Boltzmann constant, and is the temperature.
The above expression can be rearranged and written as the sum of an entropic term involving only solute coordinates and an ensemble average including also solvent enthalpy and entropy: where Δ is the configurational entropy of the solute. The above formula constitutes the basis of free energy estimation from end-point simulations detailed in the following sections. The solute vacuum potential energy was obtained by selecting the solute from the trajectories and recalculating the energy involving only solute terms.
The solute potential energy and the solvation energy were calculated from molecular dynamics snapshots of the solute. First the molecular surface was generated using the program MSMS by Sanner [38]; then the vertices and normals to the surface were read by the program Bluues [32,33] and used to compute Generalized Born radii and to compute solvation energy according to the GB model [39].
The solvent accessible surface outputted by MSMS was used to compute the apolar contribution to the solvation energy using a surface tension constant of 5 cal/(Å 2 mol) [40,41]. For entropic contributions Δ is rewritten in terms of solute conformational probability density: where (r , ) is the density in configurational space. The latter is estimated using the nearest neighbour method [42][43][44][45][46][47][48][49][50]. The basic idea of the method is to provide a description BioMed Research International 5 of the configuration in a space where a metric is defined (e.g., the torsional angle space) and then use the distance to the ℎ nearest neighbour of each sample to estimate the density of probability around that sample. When this is done with the caveats discussed in depth by Demchuk and coworkers [42], the average of the estimated density at each sample provides the estimate of the configurational entropy by the equation above.
In order to cope with the high dimensionality of the space, with largely decoupled degrees of freedom, the mutual information expansion (MIE) method [44,45] and the maximum information spanning tree method (MIST) [51,52] were proposed.
The reader is referred to the cited literature for details on the methods.
Entropic contributions due to changes in the internal degrees of freedom were calculated using the nearest neighbour and the MIST methods as implemented in our program PDB2ENTROPY (URL: https://github.com/federicofogolari/pdb2entropy) which will be described elsewhere.
Rotational and orientational entropy was computed using the nearest neighbour method in the Euclidean space approximation, which is an excellent approximation for the number of samples at hand, as described previously [50]. The calculation is implemented in our program PDB2TRENT (URL: https://github.com/federico-fogolari/pdb2trent) which will be described elsewhere. The same approach has been used by us before [53].

Molecular
Models. The X-ray structures used as starting configuration are the wild-type transthyretin (WT-TTR) taken from the RCSB Protein Data Bank (PDB) [54], PDB id: 1F41 [9], solved at the resolution of 1.3Å; the double point mutant (MT-TTR) F87M/L110M (PDB id: 1GKO [5]) resolved at 2.10Å; and the in silico engineered triple mutant F87M/L110M/S117E (3M-TTR). The structure of triple mutant was obtained from that of MT-TTR mutating the S117 by E (S117E) using the protein modeling software Swiss-PdbViewer [55]. The N-terminal residues 1-9 and C-terminals 126-127 are not present in the structure and were not modeled. Tetramers were built from deposited asymmetric units (dimers) by applying crystal symmetry operations. All the crystallization water molecules were removed prior to run the simulations.

Molecular Dynamics Simulations.
A set of independent simulation runs of different length ranging from 50 to 250 ns were performed on monomer, dimer, and tetramer of wild-type and mutant structures of TTR employing either CHARMM27 all atoms force field with CMAP correction [56] or amber99sb-ildn molecular mechanics force field [57]. The protein atoms were placed at the center of a cubic box; the system was solvated using the 3-site rigid water model TIP3P [58][59][60]. In each system an equivalent number of solvent molecules were replaced by Na + counterions to obtain a neutral system.
The systems were first minimized using the steepest descent minimization algorithm with a minimization step size of 0.1 nm and a maximum convergence force of 1000.0 kJmol −1 nm −1 . The equilibration phase was done in 2 steps: 100 ps in NVT ensemble followed by 100 ps in NPT ensemble. During the first equilibration stage, the leap-frog integrator with integration time-step of 0.002 ps was used to update the changes in the system. Particle Mesh Ewald summation [61,62] accounted for long-range electrostatics interactions. The temperature was equilibrated to a reference value of 300 K using the velocity rescaling (modified Berendsen thermostat) [63] with a coupling constant of 0.1 ps. Shortrange electrostatics and van der Waals interactions were truncated with a 10Å cutoff. All bonds were constrained with the LINCS algorithm [64]. In NPT equilibration stage, the previous parameters were still used and the pressure was stabilized to 1.0 bar using the Parrinello-Rahman pressure coupling [65,66] with a coupling constant of 2.0 ps.

Molecular Dynamics Simulation
Analysis. MD trajectories were analyzed with available structural-based tools in Gromacs-5.0.4 [41,67]. The thermodynamic stability of the tetrameric systems was further processed using the Academic License version of the bioinformatics tool FoldX [68], using the commands Stability and AnalyseComplex, respectively, to gain information on protein stability and interacting interface free energies. In all the cases, the values were averaged over the whole simulation trajectory. AnalyseComplex command outputs the Gibbs interacting free energy of binding for a complex formation (folding), say AB ( + → ), computed as ΔΔ AB = Δ AB −(Δ A +Δ B ), where Δ is the free energy of folding. Similarly the statistical effective energy function proposed by Berrera et al. [69], referred hereafter as BMF, has been used with home-written routines.
The computation of conformational entropy employs the nearest neighbour formalism whose rationale is the estimation of the local probability density around each sample by counting its number of neighbours within a hypersphere of radius equal the distance from that sample to its th nearest neighbours. Thus, the values discussed in the results section are averaged for the 10 th nearest neighbour, set as default for the routines previously cited. Conformational entropies were calculated over 5000 frames.
Pictures were either collected with PyMOL [70] or VMD [71], secondary structures were assigned using DSSP program [72], and H-bonds occupancy was computed using the readHBmap.py tool reporting only H-bonds with occupancy greater than 10 %.

Monomer to Dimer and Dimer to Tetramer Associations.
The enthalpic contributions to the association of monomers and dimers were assessed considering both single species simulations and detaching from the tetramer simulations dimers first and then monomers. In both cases favorable protein-protein interactions appear largely overestimated by the force-field used, as recently observed also for other systems [73,74]. The "enthalpic" contributions encompass the solvation energy which includes also the entropy of the solvent. For this reason the term enthalpy is used rather freely only with the purpose of separating the terms related to the Table 1: Summary of potential and solvation energy average and conformational entropy changes computed using the nearest neighbour approach from monomer (M) to dimer (D) and from dimer to tetramer (T).
WT-TTR 2M→D -4 9 . 7 ± 3.7 4.0 ± 0. potential of mean force and those computed from the analysis of conformational samples. The enthalpy of two wild-type (mutant) monomers associating in a dimer computed by the ensemble averages of the potential is -49.7 (-49.8) kcal/mol and that of two dimers associating in a tetramer is -81.0 (-77.2) kcal/mol leading to a very favorable -180.4 (-176.8) kcal/mol association enthalpy for tetramerization.
The rotational-translational entropy (Δ ) was computed for monomer-monomer and dimer-dimer association, using the Euclidean approximation for the rotational-translational space, resulting in 19.0 (17.3) e.u. for the dimer-dimer association and 17.5 (17.5) e.u. for the monomer-monomer association.
Overall, therefore, the large favorable enthalpy (likely to be overestimated here) is opposed by conformational and rotational-translational entropy which we can estimate to result in total 56.5 (61.3) kcal/mol. These results are summarized in Table 1.
The detailed analysis of the terms contributing to the enthalpy shows that covalent energy terms contribute very little to the energy difference and Lennard-Jones terms make the largest favorable contribution to association, whereas favorable electrostatic and unfavorable desolvation contributions provide an overall unfavorable but smaller (about onethird in absolute value) contribution.
The analysis provides thus a lower stability of the tetramer assembly for the mutant relative to the wild-type protein. The uncertainties of the methodology are however comparable to this difference.

Statistical Effective Energy Functions
Analysis. The structural and thermodynamic stability of both mutant structures in relation to wild-type was further probed thanks to two independent bioinformatic tools, i.e., FoldX [68] (Figure 3(a)) and the statistical effective energy function BMF [69] (Figure 3(b)). This analysis provides a complementary view of the thermodynamic picture presented above.
We computed the folding free energy differences between mutants and wild-type (ΔΔG=ΔG mutant -ΔG wild-type ) throughout the tetramer-unfolded monomers equilibrium, Figure 3.
Some words of caution are due when presenting the following data because the energy functions used take into account implicitly the entropy loss from internal degrees of freedom but not that arising from external degrees of freedom as implied by the tetramer to dimers and the dimers to monomers transitions.
Finally, notwithstanding the fluctuations and differences in underlying principles and datasets, the two approaches used here are consistent with each other as detailed below.
The panels (A) in Figures 3(a) and 3(b) point out that both mutants have positive difference in free energy as compared to WT-TTR and thus are less stable. Their average free energy ΔΔG T computed using BMF (FoldX) is +0.3 and +7.8 kcal/mol (+2.5 and +15.6 kcal/mol), respectively, for MT-TTR and 3M-TTR, confirming the decreased tetramer thermodynamic stability upon mutations, significant in the case of 3M-TTR. This result is consistent with previously determined experimental data, which indicated that monomers of 3M-TTR, at variance with MT-TTR, could not reform tetramers upon the addition of Tafamidis, a known stabilizer of the tetrameric form [24].
In order to track down all contributions, ΔΔG's values were computed and reported in subsets (B, C, and D) of Figures 3(a) and 3(b) for all the steps in the equilibrium pathway from tetramers to unfolded monomers (required before monomers self-assemble into amyloid fibrils) through dimers. We called ΔΔG D , ΔΔG M , and ΔΔG unf M the free energy for the dissociation of tetramers into dimers (T→D), for the formation of folded monomers from dimers (D→M), and for the unfolding of individual monomers (M→M unf ), respectively. Based on simulations we guessed that tetramers first break down into dimers I/II and III/IV, i.e., along the C 2 crystallographic axis, instead of I/III and II/IV. Thus, ΔΔ D = Δ D − Δ T (dissociation free energy at I/II-III/IV dimerdimer contact), ΔΔ M = Δ M − Δ D (dissociation free energy at I-II and III-IV monomer-monomer contacts), and ΔΔ unf M = Δ unf M − Δ M . ΔΔG D using BMF (FoldX) average is -1.1 and -8.7 kcal/mol (-2.5 and -26.3 kcal/mol) for MT-TTR and 3M-TTR, respectively, relative to WT-TTR. This means both equilibriums are shifted towards the right (subsets (B) in Figures 3(a) and 3(b)), i.e., the formation of dimers. The quite high (absolute) value displayed by 3M-TTR indicates the propensity of tetrameric assembly of the latter structure to dissociate into dimers I/II and III/IV as compared to the wildtype protein. Considering the dimer-monomer equilibrium (subsets (C) in Figures 3(a) and 3(b)), ΔΔG M 's average using BMF (FoldX) is -1. relative to WT-TTR. The values are however rather limited compared to the overall stability computed using BMF for the process of four WT-TTR monomers associating into 2 dimers (-14.2 kcal/mol).
Finally the change with respect to WT-TTR in folding free energy of monomers has been computed. Both BMF and FoldX predict a shift towards the formation of more stable monomers for both mutants ΔΔG unf M > 0; i.e., the  Figures 3(a) and 3(b)). The unfavorable formation of unfolded monomers in both mutants strongly correlates with the secondary structure analysis displayed in Figure 5 showing that notwithstanding the enhanced dissociation of tetramers into monomers through dimers monomers remain stable (folded). Besides, this result would strongly suggest that both monomers of our mutant variants are nonamyloidogenic, as is known as far as MT-TTR is concerned [5].

Localised Structural Transitions.
Local conformational transitions were assessed by computing the dihedral phi ( ) and psi ( ) angles of individual residues involved in the mutation and represented as Ramachandran plots, Figure 4.  Figure 4(a). These latter lie in the allowed regions for -helical secondary structure, −89 < < −39 and −66 < < −16 [75], confirming that no secondary structural transition occurred at this position.
In Figure 4(c), while both TTR variants are displaying nearly the same / averages, (-137±10, 140±9), (-136±16, 134±12), and (-131±17, 135±11), it is seen that a distinct region is being accessed by mutant variants. The region is defined by −40 < < 70 and 30 < < 140 corresponding to left handed helix conformation. The transition is observed only once for both MT-TTR and 3M-TTR within the first 10 ns of the simulation possibly as a transient adaptation to in silico mutation of residue 110.

Secondary Structural
Changes. The illustration of dynamical processes in our systems during the simulations was done by computing the secondary structure, Figure 5. The latter analysis was performed thanks to the DSSP program which determines the existence of hydrogen bonds as a criterion for the presence of secondary structure [72]. While dihedral ( , ) angles ( Figure 4) were useful for assessment of local structural changes, secondary structure is able to highlight the global structural changes. In Figure 5 we can see that the helical region exhibits large structural fluctuations and even gets disrupted from time to time, like in monomer II for all the TTR variants. F-strand undergoes some structural fluctuations at its beginning, particularly in monomers I, III, and IV. In G-H loop, some part of the structure is being converted from bend to -turn and conversely. D-strand displays quite high structural fluctuations in monomers II and IV of WT-TTR, in monomers II and III of MT-TTR, and in all monomers for 3M-TTR. A-and C-strands make the most stable parts of the molecules. E-strand in the vicinity of -helix is showing some fluctuations in monomers I, III, and IV of all the variants. In spite of structural fluctuations exhibited by major parts of the tetrameric molecular assemblies, Figure 5 clearly points out that monomers are nevertheless remaining folded; i.e., they preserve essentially their secondary structures. This observation is supported by the fact that no evident disruption and secondary structure conversion are seen in the core region (made by the -strands) of the different systems studied.

Intra-and Interchains Hydrogen Bonds Occupancy.
In order to understand the mechanism of 3M-TTR dissociation we computed the H-bonds occupancy along the MD trajectory considering both main chain-main chain, main chain-side chain, and side chain-side chain hydrogen bonds, at monomer-monomer and dimer-dimer interfaces. Figure 6 summarizes the occupancies of the interactions that were identified to significantly perturb the tetrameric assembly of TTR upon mutations. It should be noted that for 3M-TTR the occupancy reflects an average of both (starting) associated and dissociated configurations.
In Figure 6, seven major H-bonds (with ⩾ 80 % occupancy) can be seen almost symmetrically distributed in both dimers, three at I-II interface and four at III-IV interface. These are being lost in 3M-TTR with their average occupancy dropping below 40 %. These include the interactions Y114(CO)-A120(HN) and A120(HN)-Y114(O) (main chain-main chain) and T119(H 1 )-S115(O ) and S115(H 1 )-T119(O 1 ) (side chain-side chain). Most of them are located in the H-strand and some in the G-H loop.
At the interface of symmetric units (I-III and II-IV) Figure 6). Interestingly, it is worth mentioning that, in 3M-TTR, several side chain-side chain interactions are being formed following dissociation and subsequent temporary association. These involve many inner sheets residues like K15, R104, and E117, most of which show significant deviation from the starting WT-TTR X-ray structure. Analysis of hydrogen bonds connectivity confirms the pivotal role of such networks in preserving the protein integrity, typically in some dedicated portions, in this case H-strand, G-H, and A-B loops. These susceptible portions of TTR that build up the monomer-monomer and dimer-dimer interfaces are definitely playing a key role in destabilization of tetrameric subunits of TTR, especially in the case of 3M-TTR. Indeed, the disappearance of H-bonds in the previously mentioned domains along the timescale of the simulations brings a possible explanation on the integrity loss of tetramers in 3M-TTR.
3.6. Mechanism of Tetramer Dissociation. The path of dissociation of dimers in 3M-TTR may be followed during the simulation (Figure 7). The first step in conformational transition is the dissociation at the interface I/IV involving residues 17-24 and 110-123, including position 117 mutated to glutamic acid in 3M-TTR at variance with MT-TTR where the corresponding residue is a serine. The transition appears driven by the electrostatic repulsion of the pairs of acidic residues E117 close in each dimer. First the interface is weakened (up to 15 ns) and disrupted (20 ns); then both I/II and III/IV dimers remain rigid. Dimer III/IV rotates for most of the simulation about a hinge centered on salt bridges E51 (I)-R104 (III) and R104 (I)-E51 (III) and hydrogen bonds E51 (III)-T123 (I) and at the end only about the latter two interactions. The final (possibly transient) conformation is stabilized by salt bridges E117 (I)-R21 (III) and K15-E54 at I/III interface. Other interactions at II/III interface are mostly hydrophobic. No sign of dissociation or conformational rearrangement at dimers is observed during the simulation. For both dimers the proximity of E117 acidic groups should result in repulsion which is however reduced by ionic interactions.

Conclusions
Molecular dynamics simulations were used in this work to probe the structural stability and the dynamics of engineered mutants of transthyretin (TTR), i.e., the double mutant F87M/L110M (MT-TTR) and the triple mutant F87M/L110M/S117E (3M-TTR), in relation to wild-type. The analysis of trajectories reveals that mutations do not have major impact on protein structure, and the thermodynamic analysis confirms this picture.
Estimation of the free energy from end-point simulations shows that the two mutations F87M/L110M on monomers/dimer and dimers/tetramer equilibria favor the dissociation. The results are however within the uncertainties of the methodology.
Consistent with the latter analysis evaluation of stability with the statistical effective energy functions FoldX and BMF shows that the mutations and the triple mutations F87M/L110M/S117E affect differently the same equilibria and the stability of the monomers. The latter is almost not perturbed by mutations, whereas the equilibria are shifted towards the dissociated species relative to wild-type TTR. This is confirmed by the analysis of 3M-TTR which shows partial dissociation within the first 20 ns of the simulation, implying that contacts are lost at the dimer-dimer interface, whereas dimers appear fairly stable.
Overall the simulations provide a detailed view of the dynamics and thermodynamics of wild-type and mutant transthyretin and a rationale of the observed effects.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.