Probing antibody internal dynamics with fluorescence anisotropy and molecular dynamics simulations.

The solution dynamics of antibodies are critical to antibody function. We explore the internal solution dynamics of antibody molecules through the combination of time-resolved fluorescence anisotropy experiments on IgG1 with more than two microseconds of all-atom molecular dynamics (MD) simulations in explicit water, an order of magnitude more than in previous simulations. We analyze the correlated motions with a mutual information entropy quantity, and examine state transition rates in a Markov-state model, to give coarse-grained descriptors of the motions. Our MD simulations show that while there are many strongly correlated motions, antibodies are highly flexible, with F(ab) and F(c) domains constantly forming and breaking contacts, both polar and non-polar. We find that salt bridges break and reform, and not always with the same partners. While the MD simulations in explicit water give the right time scales for the motions, the simulated motions are about 3-fold faster than the experiments. Overall, the picture that emerges is that antibodies do not simply fluctuate around a single state of atomic contacts. Rather, in these large molecules, different atoms come in contact during different motions.

Introduction

We are interested in understanding the solution dynamics of antibodies. The solution dynamics of antibodies is critical to antibody function. The main body of information comes from a small number of full-length antibody structures, but static structures do not give insight into the dynamics. The internal dynamics of multidomain proteins has been a topic of interest in a number of different studies.- The focus of these studies has been on fluctuations in the solvent around the protein, the secondary structure, the domain–domain interactions and the conformations of the amino acid side chains.

Here, we study the internal dynamics of the IgG1 class of antibodies. The IgG1 class of antibodies is the most abundant human IgG subclass and the template for the majority of antibody drugs. IgG1 is composed of four polypeptide chains, two heavy chains (HC) and two light chains (LC). These four chains fold into three domains: two Fab domains that contain determinants for antigen binding and an Fc domain responsible for the effector function and binding of the Fc receptor proteins. From the early studies of antibody structure, it was apparent that Fab domains were connected to the Fc through the unstructured linker (hinge region), rendering them capable of binding to antigens separated by a range of distances., In addition, electron microscopy reports on immobilized complexes provide evidence for the significant axial rotation of the Fab domain., In subsequent years, a wide range of biophysical techniques have been used to investigate the dynamics of these versatile molecules and to look into the role of the hinge region in the structure and function of the antibodies.- Substantial efforts have been made by many investigators to obtain high-resolution X-ray crystal structures of the intact antibodies, but the lack of a stable defined solution structure of the antibodies has prevented their crystallization and only few structures of full-length antibodies are currently available., Low-resolution cryo-electron tomography studies have been very insightful in defining the conformational space explored by the antibodies.,

Previous studies suggest that the dynamics of the antibodies could play a role in their function (antigen binding, complement activation) and solution stability; however, the molecular details of these observations remain unclear.- The purpose of the present work is to look at the inter- and intra-domain motions, their timescales of motion, and to understand which amino acid residues contribute to the local and global flexibility. We do this by combining the insights obtained by time-resolved spectroscopy experiments and molecular dynamics simulations. Our MD simulations extend an order of magnitude longer than prior simulations on full-length antibodies, and we use this data to construct a Markov model characterizing the major conformational species in antibody solutions. Our results indicate that Fab and Fc fragments form multiple meta-stable protein–protein interactions, with heterogeneous protein–protein interaction surfaces based on multiple polar interactions.

Results and Discussion

Comparison of experimental and computational results

In this work, we used two independent approaches to look at the dynamics of a large multidomain biological molecule: time-resolved fluorescence anisotropy that allows determination of the rotational correlation times of the fluorescence probes site-specifically conjugated to the molecule of interest and all-atom molecular dynamics simulations that give us atomic detail structure of our system as it evolves over time in a given solution under defined conditions. To test the ability of the computational method in predicting the dynamics, we compared calculated rotational correlation times with the experimental values for three different systems of increasing size and complexity – free fluorescence probe, antibody domains and a full-length antibody.

Dye dynamics

Free dye

For a dye model, we used a common fluorescence probe Alexa 488 maleimide, which has an excited state lifetime commensurate with its rotational correlation time. This probe possesses a long (five carbon) linker between the fluorophore and a maleimide group that confers a certain degree of the flexibility to the molecule alone and once it is covalently conjugated to the protein of interest. Rotational correlation times for Alexa 488 maleimide have been reported previously and also estimated from MD simulations in two different solvents: water and methanol. We performed computational analysis of this fluorescence probe to examine its dynamics in a similar way. Our calculated anisotropy decays were best fit to the single exponential function with decay times of 71 ps and 74 ps for water and methanol (data not shown) respectively and are in good agreement with the reported values from MD simulations (51 ps and 86 ps). The observed difference between the calculated values here and those reported can be attributed to few differences in simulation parameters: (a) longer simulation times in this study (20 ns here compared with 1 ns in Schroder et al.); (b) different model for methanol used (OPLS forcefield parameters for methanol here compared with methanol parameters from GROMACS forcefield); c) all atom representation of Alexa 488 used here vs. united-atom GROMACS force field used previously. In both cases, calculated values are consistently lower than the measured parameters (170 ps and 210 ps) and the difference is due to the lower than experimentally measured solvent viscosity used in simulations. These initial results suggest that predicted rotational correlation times will also be shorter for the other two systems we examined (antibody fragments and a full-length antibody).

Dye covalently coupled to the antibody

Description of the dye behavior at the site of the conjugation to the protein is a complicated topic and a variety of mathematical models have been used to account for the movement of the probe in the measured fluorescence parameters., The reorientation time of the probe and the direction of the excitation and emission transition dipoles of the fluorophore are all parameters that can potentially have large impact on the measured values (whether these are steady-state or time-resolved values). In this work, we use three different full-length THIO-mAb derivatives of trastuzumab (humanized monoclonal antibody derivatives with engineered cysteines at defined locations; each THIO-mAb has two engineered cysteine substitutions) and their fragments (Fab and Fc) for site-specific conjugation with N-(1-pyrene) maleimide., The locations of the engineered cysteines (two different THIO-MAbs with cysteines in the Fab domain—HC121 and LC205, one THIO-mAb with cysteine in Fc domain—HC403) are shown in Figure 1. To characterize the behavior of the fluorescence probe at each of these conjugation positions, we ran all-atom MD simulations of the fragments of these THIO-MAbs with probes at each of the conjugation positions: Fab–HC121, Fab–LC205 and Fc–HC403. Starting structures for Fab trajectories were based on the crystal structure of the trastuzumab Fab (Protein Data Bank ID 1N8Z). The Fc structure was generated by removing the Fab domains from the full-length human IgG1 crystal structure (PDB ID 1HZH) up to the upper hinge lysine residue (Fc structure generated this way contained the entire intact Fc domain and all of middle and lower hinge).

Figure 1. Locations of the engineered cysteines are shown as blue spheres. Antibody domains are labeled (VL, light chain variable region; VH, heavy chain variable region; CL, light chain constant region; CH1, CH2 and CH3, heavy chain constant regions).

Each conjugation position provides a unique environment for the pyrene probe in terms of the combination of the hydrophilic and hydrophobic side chains of the protein that surround the probe. Pyrene itself is highly hydrophobic and thus expected to interact with the hydrophobic residues in its vicinity. The crystal structure of the protein provides a static picture of the environment of the probe at each conjugation position and examination of the number of hydrophobic residues around the attachment sites can serve as a qualitative measure of the extent of the dye mobility. In this work, we substituted three different residues with cysteines for site-specific conjugation: HC Ala121Cys, LC Val205Cys and HC Ser403Cys. Of these three, LC Val205 seems to be the most buried and has the highest number of hydrophobic residues within 5 Å radius, suggesting that the probe at this site would be the most rigid.

MD trajectories provide much more detailed picture of the dye environment and how it changes over the course of the simulation. Specifically, the parameters that can be estimated from the crystal structure (e.g., surface exposure, potential interaction surfaces near the dye) change over time and may not be accurately represented in the static structure. To gain a better insight into dye behavior on the protein surface, we explicitly modeled pyrene at each conjugation site. The orientation of the dye was randomly selected at the beginning of each MD run. Two Fab trajectories were collected for each of the two labeling positions on the Fab (Fab–HC121 and Fab–LC205) and one trajectory for Fc with two Fc–HC403 probes. We evaluated the mobility of the probe based on two different parameters: trajectory-averaged solvent accessible surface area (SASA) of the probe and the autocorrelation function decay time of the probe vector in the frame of the protein (see Material and Methods). Table 1 provides a summary of the trajectory-averaged SASA for all three positions. Pyrene at position HC121 seems to be the most exposed (largest SASA value) and, at position LC205, the most buried. Table 1 lists the trajectory-averaged SASA values for the conjugation residues estimated from the Fab and Fc trajectories with no dyes. All SASA values for the residues themselves are significantly lower than the corresponding values for the pyrene probe conjugated to those sites, which is a reflection of the presence of the linker between the residue and a probe. In the case of the residues themselves, Fc-HC403 is the most exposed and the other two residues display comparable surface exposure. Interestingly, two pyrenes on Fc display different surface exposure even though they are chemically identical, i.e., belong to chemically identical chains and are conjugated to the identical positions. Closer examination of the two chains in the starting structure of the Fc trajectory shows that structurally they are not identical (RMSD = 1.5 Å) and remain distinct over the course of the 100 ns MD trajectory.

Table 1. Solvent accessible surface area (nm2) of the pyrene probe or the side chain it is covalently attached to

	HC121	LC205	HC403
Pyrene-1	388	285	310
Pyrene-2	379	291	351
Side Chains	55	57	101

Side chain values were derived from the MD trajectories of the unconjugated fragments.

As expected, the mobility of the probe varies depending on the attachment position as judged by the autocorrelation function decay time in protein frame (Table 2). Pyrene at HC121 position appears to be the most mobile in our simulations. As in the case of surface exposure, two pyrenes on Fc domain (position HC403) remain distinct in terms of reorientation times, with one being faster than the other (68 ns vs. 1028 ns). In each MD run, we observe pyrene probe forming transient interactions with the hydrophobic patches in the vicinity of the attachment site. In the case of the two Fc probes, different starting orientations were chosen for the MD run, as reflected in the surface exposure (Table 1). As a result, HC403-1 forms a long-lived interaction with antibody side chains resulting in much longer rotational correlation time in the protein frame.

Table 2. Computationally derived pyrene dye rotational correlation times predicted from MD data.

	HC121	LC205	HC403
Pyrene-1	12 ns	423 ns	1017 ns
Pyrene-2	7 ns	127 ns	68 ns

HC121 and LC205 data are from two different trajectories. HC403 data are from a single trajectory with two probes (protein frame).

Antibody dynamics

Fab and Fc dynamics

Modeling pyrene probes in MD trajectories allows direct comparison of the hydrodynamic properties of the antibody fragments as measured experimentally with the predicted parameters. Fab simulation was performed in a 850 nm3 box with approximately 25,500 water molecules and in the case of Fc − 1,000 nm3 box with approximately 33,700 solvent molecules. Simulations were performed at pH 7.0 and 25°C. At this pH, the trastuzumab Fab is positively charged (+5). Chloride ions have been added to the simulation box to make the system neutral. The Fc is neutral at this pH.

Experimental and computational results for rotational correlation time of the antibody fragments are summarized in Table 3. In all cases, experimental anisotropy decay data was fit to the sum of two exponentials and an average anisotropy decay time was calculated (see Material and Methods). As can be seen from Table 3, the experimental results give approximately the same rotational correlation time for both conjugated positions on the Fab, whereas the Fc gives a slightly faster rotational correlation time. As the Fab and Fc volumes are comparable, we attribute the differences between Fab and Fc rotational correlation rates to local flexibility of the protein around the attachment position or an altered intrinsic flexibility of the fluorescence probe at each position (as defined by distinct surface exposure and dye-side chain contacts). In addition, the extent of fluorescence depolarization depends on the orientation of the emission transition moment of the fluorescence probe relative to the axis of rotation of the macromolecule to which it is conjugated. Previous reports on time-resolved anisotropy of the macromolecules attempt to extract individual contributions of different sources of depolarization (e.g., dye, domain and whole protein tumbling); however, because the independent determination of individual contributions is not possible with the experimental tools that we have and the problem of defining these contributions from a single decay does not have a unique, well-defined solution, use of the single (or the average) decay time is an appropriate measure of dynamics for a given labeling position.

Table 3. Comparison of experimental and molecular dynamics.

	Fluorophore Conjugation Position: Fab and Fc
Method	Fab–HC121	Fab–LC205	Fc–HC403
Experimental	20 ± 3 ns	24 ± 2 ns	11 ± 3 ns
MD (lab frame)	3.4 ± 0.4 ns	11.1 ± 2.3 ns	11.2 ± 3.8 ns

Rotational correlation times for trastuzumab fragments (Fab and Fc) determined using time-resolved anisotropy decay or calculated from the MD trajectories (lab frame).

As can be seen from Table 3, predicted rotational correlation times from the MD data are in overall agreement with the experimental measurements, but there are obvious differences,. As mentioned before for the free dye, viscosity of the water model used here is lower than the experimental viscosity, leading the modeled dynamics to be faster than the observed values (with the exception of HC403). In addition, we hypothesize that the local protein environment of the probe (based in this case on the crystal structure) may not necessarily be the most representative state found in freely diffusing molecules. Rather, it is likely (for at least some of the dye conjugation positions) that local conformation determined by crystallography or a given local structure resides in a local energy minima, thus biasing the dynamics of the probe for that specific local environment. As a result, we observe that the apparent tumbling time of the Fab—as predicted from pyrene at position HC121—is faster than the other two positions, a result of the greater flexibility of the dye itself at this position. To gain a better understanding of the dye dynamics, starting structures based on NMR data would probably offer an advantage; however, no such structures are available for the antibodies due to their large size.

Rotational correlation time can also be predicted from the autocorrelation function decay of the bonds of the dye attachment residue; specifically, Cα–Cβ and Cα–N bonds of these residues (Table 4). These values were determined from the trajectories of the fragments without probes and trajectories with probes present. In both cases, rotational correlation times are similar (the presence of the dye does not appear to affect the local dynamics of the bonds) and compare well with the experimentally measured parameters (compare with the experimental values in Table 3). In this case, the determined rotational correlation time does not explicitly incorporate dye dynamics at each specific position, but it does offer a good approximation of the local dynamics as judged by the similarity of the calculated and experimental values. From the previous points, it follows that measuring dynamics from explicitly modeled probes and from the bonds of the residues gives similar results and either can be used to predict the dynamics.

Table 4. Rotational correlation times calculated for the bonds of the given resides from the MD data of the antibody fragments.

	Fab–HC121		Fab–LC205		Fc–HC403
Trajectory	Cα–Cβ	Cα–N	Cα–Cβ	Cα–N	Cα–Cβ	Cα–N
with dyes	13 ± 2 ns	16 ± 2 ns	16 ± 2 ns	12 ± 3 ns	15 ± 2 ns	16 ± 1 ns
no dyes	10 ± 2 ns	16 ± 0.1 ns	13 ± 4 ns	11 ± 1 ns	16 ± 7 ns	15 ± 1 ns

Trajectories with or without fluorescence probe present were used.

Dynamics of full-length antibody; all atom molecular dynamics simulations

To test the ability of molecular dynamics one step further, we compared the parameters obtained from 1 μsec (total) of all-atom MD simulations of full-length antibody with experimentally measured parameters.

Full-length trastuzumab is computationally challenging to model with molecular dynamics simulations because of its large size (1328 residues, in contrast with commonly modeled proteins that have ~50 residues or less) and multiple domains (almost all modeled proteins are single domain, whereas trastuzumab has three distinct domains connected via flexible linker). A critical challenge is to obtain enough sampling of the conformational dynamics. We address this issue by comparing our computational results with fluorescence anisotropy.

We maximize conformational sampling in our simulations by the following two approaches. First, we run our simulations on a multi-processor Linux cluster using optimized molecular dynamics software (multiprocessor GROMACS 4 with virtual hydrogens, allowing for 4 fs timesteps (see Materials and Methods for more details). We collected approximately 2 μs of data, an order of magnitude more than has been previously published for full-length antibodies. Second, we use multiple starting structures obtained by (a) homology models from two different templates (see below), and (b) coarse-grained modeling with an anisotropic network model to calculate the extremes of the lowest frequency normal modes (see Materials and Methods for details). Simulations where performed in explicit water in a 43,000 nm3 box. The box contained 6,870 protein atoms and approximately 135,000 molecules of solvent.

Initial full-length trastuzumab homology models were generated based on two crystal structures available in PDB public database, human IgG1 (PDB ID 1HZH) and mouse IgG1 (PDB ID 1IGY) (Fig. 2). The templates are structurally distinct from one another in the following two ways: (1) the Fab fragments are rotated ~180° along its long axis in one structure vs. the other, and (2) in one template there is significant inter-domain interaction between one Fab and the Fc, whereas in the other there is less interaction between the fragments and the Fc is perpendicular to the plane of the 2 Fabs. The first model was based on template 1HZH, a IgG1κ human antibody from an HIV immune patient. The second model was based on template 1IGY, a mouse anti phenobarbital antibody. 1N8Z trastuzumab bound to the extracellular domain of the HER2 receptor was used as the Fab template for both homology models. The G2 form of sugar was used, with coordinates taken from 1L6X, the crystal structure of rituximab Fc. The two structures differ significantly so we chose to create both structural models because there is no compelling evidence for the prevalence of either one in solution. Anisotropic network model (ANM, see Materials and Methods) was used to generate starting structures for a total of 12 MD trajectories (six for 1HZH-based models and six for 1IGY-based models) 70 nsec to 500 nsec in length. All starting structures are pictured in . summarizes all the simulations we ran. In the case of full-length antibody simulations, we did not model fluorescence probes. The parameters describing protein dynamics were calculated from the bonds of the corresponding residues.

Figure 2. Crystal structures of human (A) and mouse (B) IgG1 used for generating homology models of trastuzumab . In both structures, Fab fragments are structurally distinct and have been treated as distinct domains throughout the manuscript.

Measured anisotropy of full-length trastuzumab

Results for the measured rotational correlation times are summarized in Table 5. For the full-length mAb, each conjugation position gives a different value of the rotational correlation time. We attribute these differences primarily to the unique environment of the probe. The Fab–LC205 shows a slower rotational correlation time as compared with Fab–HC121, reflecting a more constrained environment in the full-length mAb. This difference was not apparent in the measured anisotropy of the fragments, demonstrating that the local environment between the two positions differs between the fragment and a full-length antibody.

Table 5. Experimental (first three rows) and calculated (MD) rotational correlation times for the full-length trastuzumab as determined from three different conjugation positions

	Fluorophore conjugation position: Full length mAb
Solution conditions	Fab–HC121	Fab–LC205	Fc–HC403
pH 7.0, no salt	46 ± 4 ns	73 ± 7 ns	35 ± 6 ns
pH 7.0, 150 mM NaCl	44 ± 5 ns	73 ± 5 ns	33 ± 4 ns
pH 7.0, 450 mM NaCl	49 ± 4 ns	75 ± 5 ns	37 ± 5 ns
MD (lab frame)	41 ± 24 ns	43 ± 21 ns	47 ± 16 ns

To better understand these differences, we also performed anisotropy measurements under three different salt concentrations to assess the impact of ions on the dynamics of the antibody. As can be seen from Table 5, salt concentration does not have an effect on the apparent rotational correlation times, and thus the apparent hydrodynamic radius of the molecule remains unchanged. In addition, as a measure of sensitivity of time-resolved anisotropy to the changes in the dynamics of the macromolecule, we determined the rotational correlation times for the trastuzumab THIO-mAb that lacks hinge disulfide bonds (to achieve this, hinge cysteine residues have been substituted to serines; we refer to this construct as “hingeless”). The fluorescence probe in this case was located on the Fab fragment (LC205). Absence of the covalent link between the two heavy chains should allow more flexibility to the Fab fragments and effectively decouple the Fab motion from the rest of the antibody. As expected, for this molecule, we see the decrease in the correlation time (see Table 6). The value of the rotational correlation time for the Fab–LC025 within the “hingeless” constructs is intermediate between the free Fab and a Fab in full length trastuzumab and thus it appears that Fab motion does not become completely decoupled from the rest of the antibody. The same table shows the values of rotational correlation times for antibodies derived from two different sources (E.coli and CHO cells). The only difference between the two is their glycosylation state (E.coli-derived antibodies lack sugars). The rotational correlation times for these molecules are the same, suggesting that the absence of the glycosylation does not substantially alter the apparent hydrodynamic radius of the molecule. Reduction of the intact antibody (under these conditions only interchain disulfides are reduced) results in the value of the rotational correlation time comparable to the one for the “hingeless” derivative (second row, Table 6).

Table 6. Rotational correlation times for trastuzumab LC205 full lenth antibody derived from E.coli or CHO cells (intact antibodies) and trastuzumab derivative with hinge disulfides removed (“hingeless”), in native and reduced states

	Full length LC205
Solution conditions	E.coli	E.coli “hingeless”	CHO
pH 7.0, 150 mM NaCl	73 ± 5 ns	54 ± 6 ns	74 ± 9 ns
DTT, 12 h	53 ± 7 ns	49 ± 6 ns	50 ± 5 ns

Comparison of measured anisotropy vs. calculated

The last row of the Table 5 shows the calculated rotational correlation times based on the MD data for the full-length antibody. These values represent an average of all 12 trajectories based on the two initial models since the decay times calculated from trajectories based on these models were comparable (data not shown). It appears that the dynamics as observed using MD simulations is not affected by the choice of the initial structural model.

As mentioned previously, the calculated rotational correlation times were estimated from the decays of the autocorrelation functions of the Cα–Cβ and Cα–N bonds of the corresponding residues (and not explicitly modeled probes). As a result, the predicted values lack the contribution from dye dynamics. Despite this, calculated values are in excellent agreement with the measured correlation times, again justifying the use of the computational tools in describing solution dynamics of large multi-domain proteins.

Computational dynamics of trastuzumab: Correlated motions

Having validated our computational models with fluorescence anisotropy, we next describe the conformational dynamics of trastuzumab since correlated motion in proteins is essential to function. To quantify the correlations, we use mutual information (MI), a metric from information theory that captures all correlations including nonlinear and higher order correlations. We use MI as a metric because MI in molecular dynamics has been shown to capture biologically relevant dynamics,,- and a comparative analysis revealed that correlations captured by MI were more functionally relevant than methods such as principal component analysis.

MI between all residue pairs was calculated and normalized as described in Materials and Methods. Briefly, the trajectory of each residue was represented by the distance of the backbone Cα atom from its moving average position (the average position in a 500 ps window centered on the current time step). MI is normalized by the joint entropy H(X,Y) so that MI between lower entropy residues is weighted equally with MI between higher entropy residues. In this report, normalized MI, MInorm, are provided. Figure 3 shows the MInorm between all 1328 × 1328 residue pairs in the following order: Fab fragment 1 (light chain 1, and the variable and CH1 domain from heavy chain 1), Fab fragment 2 (light chain 2, and the variable and CH1 domain from heavy chain 2), and the Fc fragment (CH2 and CH3 from heavy chain 1, and CH2 and CH3 from heavy chain 2). The rows and columns of the MI matrix in Figure 3 are the residues, and the color-coded values in the matrix correspond to the amount of correlation seen in the dynamics as quantified by MInorm. Strong correlated motion between two residues (high MInorm) is red, and two residues that move independently (low MInorm) are blue. The definition of MI results in the following two properties of the matrix: (a) the highest MI is between a residue and itself [MI(X,X)], the matrix diagonal (perfect correlation), and (b) the matrix is symmetric because MI(X,Y) = MI(Y,X). In all figures, we omit the values on the diagonal to redistribute the color scale for clarity (instead of the scale going from 0 to 1, the scale goes from 0 to the maximum non-diagonal MInorm value). Strong correlated motion in residues that are adjacent in primary sequence is seen as the thickness of the diagonal. The secondary structure of trastuzumab is largely antiparallel β strands arranged in β sheets, visible as high MI lines perpendicular to the diagonal because residue X in one β strand has high MI to residue Y in an adjacent β strand, residue X + 1 has high MI to residue Y − 1, etc. As expected, there is a strong correlation in motion between residues that are adjacent or interacting.

Figure 3. Correlated motion between all residue pairs in trastuzumab is represented as a ‘heat map’. In this ‘heat map’, the x- and y-axes identify the residue numbers and the color indicates the amount of correlation motion in the molecular dynamics simulations (red, high; blue, low). Correlation is calculated by MInorm (a general form of correlation that takes both linear and nonlinear correlation into account). MInorm ranges from 1 (maximum correlation; e.g., MInorm between a residue and itself) and 0 (two residues that move independently). To expand the color scale to clarify the figures, the diagonal line representing MInorm between a residue in itself is not included in the color scale. The Fab and Fc fragments are easily identifiable as four Ig domains in which intra-Ig domain residues have highly correlated dynamics.

Correlated motion reveals tight coupling within the Ig domains

MI reveals a high degree of coupling within each Ig domain (Fig. 4; each of the 12 Ig domains are visible as bright squares in the matrix). Residues within each of the 12 Ig domains are most coupled to other residues within the same Ig domain and there is tight coupling throughout most residues in the domain. Secondary structure interactions between the β strands of the β sheets are the most tightly coupled (e.g., interactions between anti-parallel strands are visible as lines perpendicular to the diagonal). This suggests that residues within an Ig domain move as a unit.

Figure 4. Average correlated motion between the 12 Ig domains. (A) Residues within an Ig domain move in a correlated fashion (Fig. 1), and we average all intra-Ig domain residues to highlight correlated motions between Ig domains. Each square on the diagonal is the average correlated motion between residue pairs in the same Ig domain. Off-diagonal squares indicate correlated motion between Ig domains. High correlated motion (high MInorm) is red; low correlated motion (low MInorm) is blue. (B) Cartoon representation of the correlated motion between Ig domains. The highest correlation is seen between intra-Fab variable Ig domains and between CH3 Ig domains. Intra-Fab CH1 Ig domains have correlated motion that is slightly lower than the variable and CH3 Ig domains. CH2 domains, which are glycosylated, move relatively independently (low correlated motion). With the exception of CH2, Ig domains in the “horizonal” direction (protein–protein interaction along the direction of the β strands) move in a correlated fashion and Ig domains in the “vertical” direction (domain–domain interactions between the loops at the ends of β strands) do not move in a correlated fashion. Correlated motion between the two intra-Fab Ig domains could facilitate binding to antigen, and the lack of correlated motion between the variable Ig domains and CH1 insulates antigen binding from effector functions.

Whereas there is no direct experimental evidence available for coupled motion in Ig domains, the following properties of antibodies are consistent with tightly coupled dynamics between residues in Ig domains. Antibodies are highly thermostable. Ig domains are modular, which can be seen in the following two properties. First, many immune proteins are composed of multiple Ig domains in various configurations. Second, proteins composed of subsets of the Ig domains of full-length antibodies (e.g., Fab alone, or just the variable Ig domains that contain the antigen binding sites) maintain similar antigen binding function. Finally, Ig domains are structurally conserved across proteins,, as well as between various states (e.g., the RMSD between trastuzumab bound to antigen or alone is 2.5 Å).

Hinge residues move independently and have the greatest phi-psi local conformational entropy

The hinge region acts as a flexible linker between the fragments. Residues in the hinge have been grouped into three regions according to their relative orientation to the inter-chain disulfide bonds: the upper hinge consists of residues between the disulfide bonds and the Fab fragment, the middle hinge consists of residues near the disulfide bonds, and the lower hinge consists of residues between the disulfide bonds and the Fc fragment.

MI of residues in the hinge region reveals that these residues move independently from the rest of the antibody (Fig. 5A), consistent with their characterization as a flexible linker between fragments. To further characterize flexibility, the phi-psi local conformational entropy of every residue was calculated (see Materials and Methods for more details). Briefly, the Shannon entropy in phi/psi angles was calculated for all residues of the full-length antibody. Residues in the hinge region were the most flexible (exhibited the highest phi/psi entropy), along with residues at the N- and C-termini. Figure 5B shows the conformational entropy of all residues in the hinge region. Residues with high conformational entropy are colored red, and residues with lower conformational entropy are colored yellow. Entropy clearly reveals the three regions of hinge residues—high entropy upper hinge, low entropy middle followed by high entropy lower hinge. The borders of these three regions are mostly created by the proline residues (this type of “mosaic” structure has been observed in NMR experiments).

Figure 5. The hinge region in trastuzumab is the most flexible region in the molecular dynamics simulations. Differences in dynamics delineate an upper, middle, and lower hinge region. (A) Correlated motion between all residue pairs in the trastuzumab hinge region represented as a ‘heat map’. Residues with high correlated motion are colored red (high MInorm); residues that move independently are colored blue (low MInorm). Motion in hinge residues is less correlated with adjacent residues than Ig residues (‘thinner’ diagonal), and hinge residues move relatively independently from the Ig domains (low MInorm between hinge and Ig residues). The middle hinge region, in which residues in one heavy chain are disulfide bonded to residues in the other heavy chain, has higher MInorm between adjacent residues relative to the upper and lower hinge (middle hinge residues have a ‘thicker/slightly orange’ diagonal relative to upper and lower hinge residues). (B) Residues in the hinge region have the highest conformational entropy. Conformational entropy is indicated by color on residues in the hinge. Red indicates high conformational entropy (high entropy in phi/psi angles), and yellow indicates medium conformational entropy. Residues in the upper and lower hinge region have the highest conformational entropy.

No common dynamics in complementarity determining regions

Each Fab fragment binds to an antigen via six loops called complementarity-determining regions (CDRs). MI of heavy chain CDRs and conformational entropy of the heavy and light chain CDRs are shown in Figure 6. Each CDR loop exhibits different dynamics. Heavy chain CDR3 residues have the lowest MI between them and the highest conformational entropy. This indicates that the residues in this loop move relatively independently from one another and with high flexibility. The residues in this loop also move independently from the rest of the antibody (low MI between residues in the loop and residues in the rest of the antibody). The most tightly coupled loops (high MI between residues in the loop, indicating that the residues in the loop move together) are light chain CDR1, light chain CDR2, and heavy chain CDR1, consistent with the canonical structure., Some loops have motion that is coupled to non-CDR residues in the Fab domain (high MI between one or two residues in the loop and the other Fab residues). The highest MI interactions between heavy chain CDR loops and non-CDR regions are shown in Figure 6. We calculated the phi-psi conformational entropy of the Fab domain to further characterize flexibility, as we did for the hinge region (Fig. 6B). The CDR loops have the greatest conformational entropy within the Fab domain, consistent with their previous characterization as structurally flexible. Heavy chain CDR-3 has the greatest conformational entropy of the six CDR loops.

Figure 6. Heavy chain CDR-3 moves independently and has the highest conformational entropy, consistent with previous data showing that this CDR shows the greatest variability. (A) Correlated motion between all residues in the trastuzumab heavy chain variable Ig domain. Residues with high correlated motion are colored red (high MInorm); residues that move independently are colored blue (low MInorm). The three heavy chain CDR regions are circled. Heavy chain CDR-3 moves the most independently of all CDR regions (low correlated motion, light chain variable Ig domain not shown). (B) Conformational entropy (phi/psi entropy) in the variable Ig domains indicates that the greatest entropy is seen in the CDR regions. Heavy chain CDR-3 has the greatest conformational entropy of all CDR. Red indicates high conformational entropy (high phi/psi entropy), and blue indicates low conformational entropy (low phi/psi entropy). regions.

Correlated motion reveals interactions between Ig domains

Having established that residues within an Ig domain move in a coupled fashion, we reduce the 1328 × 1328 residue MInorm matrix to a 12 × 12 matrix where each element in the matrix represents the average MInorm between all residues in the corresponding Ig domain(s) (Fig. 4A). The MI between Ig domains are more easily viewed in this reduced matrix.

Within Fabs/Fc in general, there is coupled motion between Ig domains that interface “horizontally” along their β strands and less coupled motion between Ig domains that interface “vertically” between their loop regions (Fig. 4B). For example, the two variable domains within the Fab fragments move together (interchain), and move more independently from their adjacent CH1 domains (intrachain). The exceptions to this are the CH2 domains, which do not move in a coupled fashion. These domains are glycosylated, and the sugars are the interaction site between the domains, resulting in reduced coupling between the domains. Figure 4B summarizes the coupling between the Ig domains.

As in the above discussion, there is no direct experimental evidence for coupling between Ig domains; however, existing data are consistent with the coupling patterns in the computational data. In the case of more coupling between the two variable domains than between the variable domains and their adjacent CH1 domains, this data can be viewed in the context of the antigen binding function of the variable domains. Antigen binding occurs across three loops (CDRs) per variable domain for a total of six CDRs per Fab fragment (not all loops may be involved in binding). Thus, binding is coordinated across two domains, which would be facilitated by coupled dynamics between them. In addition, the variable domains move relatively independently from the rest of the antibody, which may decouple the antigen binding function from the effector function of the Fc fragment. The lack of the coupled motion between the glycosylated CH2 domains is consistent with the fact that their interaction is primarily mediated through the two sugars attached to Asn297 and not between the residues of the domains. Carbohydrate-carbohydrate association is known to be weak relative to domain–domain interactions, and CH2 domains have higher crystallographic b-factor values than the adjacent CH3 domains, suggesting greater flexibility. The CH2 domain is known to melt at a lower temperature, suggesting the lack of interdomain residue interactions that may result in reduced stability of the CH2 domains relative to other domains.

Between Fabs/Fc, there is coupled motion between the “inner” Ig domains (the CH1 and CH2 domains), whereas the “outer” Ig domains (variable and CH3 domains) move relatively independently from the rest of the antibody (Fig. 4B). We investigate the molecular details of these interactions in the next section.

Interactions between fragments are largely polar and non-specific

To investigate the nature of the coupled motion between fragments, we examined each trajectory in more detail. In every trajectory, there is at least one occurrence of inter-fragment Ig domains in contact (contact defined as within a 6 Å distance). Both Fab–Fab and Fab–Fc contacts are made exclusively through the “inner” Ig domains (CH1 and CH2). Examination of the interaction surfaces in detail revealed that the interaction surfaces are heterogeneous (and thus interactions are non-specific), with largely polar interactions mediated through multiple salt bridges. For example, the 500 ns trajectory started from the 1HZH homology model (trajectory 0), and a 90 ns trajectory started from the same homology model (trajectory 2); both have interactions between the two Fab fragments (between the CH1 domains). The interaction interfaces in these two trajectories are different. In trajectory 0, the interaction occurs via the following salt bridges: ASP336 and LYS188, GLU187 and ARG425, GLU187 and LYS404, GLU213 and LYS404, GLU401 and LYS183, GLU427 and ARG211, and GLU427 and LYS190 (Fig. 7). In trajectory 2, the interaction occurs via the following salt bridges: ASP151 and ARG425, ASP336 and LYS149, and GLU195 and LYS340 (Fig. 7). These interactions are largely through loop regions, but salt bridges to secondary structures also occur. These heterogeneous interactions sometimes persist over the course of a given trajectory and sometimes are transient.

Figure 7. Salt bridges between fragments form and break during the simulations, and the domain–domain interfaces are heterogeneous (different trajectories form different salt bridges, even when the same fragments are interacting). A cartoon representation close up of the two Fab fragments indicates residues that participate in inter-Fab salt bridges for trajectory 0 and trajectory 2 (shown as a representative example). Residues that participate in salt bridges in trajectory 0 only are shown as cyan spheres. Residues that participate in salt bridges in trajectory 2 only are shown as purple spheres. Residues that participate in salt bridges in both trajectories are shown as blue spheres.

All simulations discussed to this point were performed with no salt present and resulted in polar contacts forming between different surfaces of the antibody. To test the validity of these observations for physiological salt concentrations, we also performed MD runs in the presence of 150 mM NaCl using starting models based on both 1HZH and 1IGY structures. We collected a total of six trajectories, each 160 ns long. Again, as in the case of no salt simulations, we observe that configurations, having at least one salt bridge, between any two of the domains, exist. We also observe contacts in which no salt bridges between any of the domains exist. In an attempt to understand inter-domain (Fab–Fab/Fab–Fc) salt bridges of trastuzumab in solution, we compared the no salt simulations to the salt simulations by calculating the average number of salt bridges between each of the domains and determining if salt has an effect on the average number.

Table 7 lists the average numbers of inter-domain salt bridges. The numbers of salt bridges are larger for “salt” simulation averages than for the “no salt” simulation ones. The differences are not significant except for the Fab1–Fab2 inter-domain salt bridges from the simulations with starting models built from 1HZH template. In theory, adding salt will weaken inter-residue electrostatic interactions. In this case however, it appears that salt screens the electrostatic repulsion between the Fab fragments (each Fab has a net +5 charge at pH 7.0) and does not significantly affect the Fab–Fc interaction (Fc is not charged at this pH). We thus conclude that domain–domain interactions (with or without polar contacts) exist under both conditions tested. provides a summary of the total number of contacts formed over the course of salt and no salt trajectories.

Table 7. Average number of salt bridges between trastuzumab fragments in the context of the full-length antibody

Starting structure	Interacting fragments	Number of salt bridgesat 3.2 Å		Number of salt bridgesat 5Å
				No salt	Salt	No salt	Salt
1HZH	Fab1–Fab2	1.417	3.051	2.934	6.89
	Fab1–Fc	0.286	0	0.646	0
	Fab2–Fc	0.675	0.936	1.159	1.244
1IGY	Fab1–Fab2	0.549	0.562	1.135	1.093
	Fab1–Fc	0.518	0.61	1.641	1.633
	Fab2–Fc	0.45	0.476	0.83	0.802

Trajectories with initial structural models based on both crystal structures (human and mouse IgG1) were used.

Experimental evidence for inter-fragment interaction comes from structural studies of antibodies in crystallized (X-ray crystallography) or immobilized (electron microscopy) form. X-ray crystallography of a full-length human IgG antibody (used as a template to build a homology model for this study) reveals significant interaction between one of the Fab fragments and the Fc fragment. X-ray crystallography of multiple proteins with IgG domains reveals a diversity of interaction surfaces between two Ig domains or between Ig domains and other proteins, and indeed each part of the surface has been part of a domain–domain interaction interface. Electron microscopy also reveals conformations that are consistent with inter-fragment interactions,, albeit in less detail than X-ray crystallography. Antibodies are soluble proteins, and, as such, their surface is polar, which is consistent with polar interactions between fragments. Finally, there is a high effective local concentration of the Fab/Fc fragments in the vicinity of each other. This high concentration, estimated to be millimolar, likely drives these non-specific interactions. Finally, the flexible hinge region connecting the fragments allows for multiple relative orientations of the fragments exposing different surfaces for potential non-specific inter-fragment interaction.

The solution state of trastuzumab is a heterogeneous population of conformations with inter-fragment interactions

To estimate contact rates between trastuzumab domains from the simulation data, we use a Markov State Model (MSM) approach. An MSM is a kinetic model that describes conformational dynamics as transitions between kinetically metastable states. If such a set of suitably Markovian states can be identified, estimates of transition rates between states can be used to obtain a complete description of the thermodynamics and kinetics of the system.- MSM approaches have been very useful, for example, in integrating large-scale simulation data to build networks of protein folding at long time scales.,

Here, we used the MSM framework to address the more modest task of estimating the timescales of motion involved in interdomain contact formation. Our goal is to obtain simple estimates of the conformational dynamics of trastuzumab in solution, rather than detailed conclusions about specific mechanisms of domain interaction or specific kinetic rates. To do this, we first manually define a set of eight conformational states, each having a unique combination of interdomain contacts. We then estimate rates of making/breaking interdomain contacts from the simulation data (taking special care in estimating statistical uncertainties from finite sampling effects). The resulting 8 × 8 rate matrix is diagonalized to obtain a spectrum of implied timescales of motion, as described in more detail below.

We manually chose a state decomposition of 23 = 8 states, one for each unique configuration of possible inter-domain contacts (Fab1–Fab2, Fab1–Fc, and Fab2–Fc; Fig. 8A). Trajectory snapshots were assigned to states based on the number of atomic contacts between domains (two atoms are considered in contact if they are within 6 Å of each other). A threshold between 500 and 800 atomic contacts is used to determine state transitions: if two states are in contact, they must have less than 500 atomic contacts before they are considered to lose contact, and if two states are not in contact, they must have more than 800 atomic contacts before they are considered in contact. This procedure serves to prevent the overestimation of transition rates due to multiple barrier crossings. Although the states were chosen manually, the decomposition was chosen with the goal of having minimal assumptions about the molecular mechanisms underlying contact formation. While inter-domain contacts may form at a number of different interaction sites, it is reasonable in this case to consider a simple coarse-grained model that assumes an average contact rate for each pair of domains, integrated over all possible interactions.

Figure 8. (A) A schematic diagram of the three domains of trastuzumab. (B) A visual depiction of the transition matrix (T) (for a lag time of τ = 1 ns) compiled from observed transition counts across all 12 simulation trajectories. (C) Transition probabilities from T, shown as a graph. Transition probabilities greater than 0.02 are shown as a solid arrow.

In keeping with this, we analyzed 12 independent simulations started from different initial conformations. The starting configurations came from homology models built from immunoglobulin crystal coordinates (PDB IDs: 1IGY, 1HZH) and coarse-grained ANM models. Because of the large size of trastuzumab, the simulations do not ergodically sample all possible contact states. In many cases, different simulation trajectories sample disjoint regions of state space, due to their very different starting conformations. Thus, there is a great amount of sampling bias that enters into the MSMs we construct. A nonparametric bootstrap analysis is used to estimate the variance in contact rates introduced by this finite sampling bias (see Materials and Methods).

Figure 8B shows the transition matrix T (for a lag time of τ = 1 ns) compiled from observed transition counts across all 12 simulation trajectories. The transition rates show that most of interconversion between states is predicted to arise from the making and breaking of Fab1–Fc and Fab2–Fc contacts. The equilibrium state probabilities (Fig. 8C) predicted by the MSM show that states with one or two contacts are the most likely.

The eigenvector structure of the slowest-timescale mode shows that the slowest contact events correspond to the making and breaking of Fab1–Fab2 contacts (Fig. 8C). The posterior distribution of the slowest implied timescales calculated by our bootstrap procedure predicts a mean implied timescale of 35.5 ns (Table 8), although the range distribution is quite broad, ranging from 7.1 ns to ~2.0 µs (the smallest and largest values sampled in the bootstrap, respectively).

Table 8. Mean implied timescales for each relaxation mode in the MSM (for a transition matrix T compiled from observed transition counts across all 12 simulation trajectories)

Mode	Mean implied timescale (ns)	Standard deviation (ns)
1	35.548	32.240
2	11.252	3.612
3	7.897	2.490
4	2.818	1.384
5	0.489	0.689
6	0.154	0.036
7	0.141	0.023

Based on these results, we conclude that in solution, trastuzumab has a heterogeneous population of conformations, most with one or more inter-fragment interactions. Although accurate estimates of dynamical timescales are limited by finite sampling, our best estimate from the simulation data are that these events occur on the ~40 ns timescale, although motions could be as slow as microseconds. Hence, we expect trastuzumab in solution to be interconverting between these metastable states.

Materials and Methods

Molecular dynamics simulations

Molecular dynamics (MD) simulations of the complete trastuzumab antibody and its Fab fragment were performed with the GROMACS 4.0 simulation software package. The simulations were performed on a 256 processor computer cluster running the Linux operating system with Infiniband® network connectivity.

Homology models were built and used for the initial structures of the complete trastuzumab antibody. The templates for the homology model were the human IgG1 (1HZH) and murine IgG1 (1IGY) protein crystal structures. Homology models were built by structurally aligning the Fab crystal structure (1NZ8) to the respective template Fabs and performing a direct coordinate replacement. Missing atoms, in the template structure, were built in using the SWISS-MODEL- homology modeling server. In particular, residues in the upper hinge of one heavy chain were missing in the template pdb 1HZH. In addition to missing residues, the hinge region of template 1HZH was missing one of two interchain disulfide bonds. The hinge coordinates were relaxed to within disulfide bond distance using energy minimization. This reconstructed hinge was used to replace the hinge region for template 1IGY, a murine antibody with a different hinge.

The structure of the two trastuzumab homology models (see e.g., Fig. 2) were further analyzed using an anisotropic network model (ANM), to calculate the normal modes of the trastuzumab antibody (http://ignmtest.ccbb.pitt.edu/cgi-bin/anm/anm1.cgi). A 15-Å cutoff was used to define interacting residues. Eight structures, four for each of the homology models, along the two lowest amplitude modes were selected as additional starting structures for the MD, for a total of ten initial starting structures. Structures that represented the extremes of the range of motion were chosen as starting structures to maximize exploration of conformational space.

All trastuzumab structures were glycosylated in G2 form with sugars on HC ASN300 in the Fc domain.

The 1N8Z trastuzumab Fab crystal structure was used as a starting structure for Fab only simulation. The Fab was conjugated with N-(1-pyrene) maleimide, which forms a succinimide upon conjugation, at LC resides 205 and HC residue 121 in some of the Fab only simulations. The initial starting configuration was by hand followed by equilibration and randomization.

The OPLSAA force field, was used to model the protein. The charge state of the titrateable residues was evaluated using the empirical method PROPKA., All the residues were set to their canonical protonation state. The OPLS carbohydrate force field, was used to model the bonded and Lennard–Jones parameters of the carbohydrate used to glycosylate the full-length antibody. The charges on the carbohydrate were obtained using the semi-empirical charge model AM1-BCC as implemented in the Antechamber software package., The boned and Lennard–Jones parameters of the N-(1-pyrene)maleimide were obtained using the OPLS forcefield where available. The charges on the N-(1-pyrene)maleimide were obtained AM1-BCC.

OPLS atomstypes were used for the Alexa 488 atoms. The charges were calculated using antechamber. Simulations were performed as described above. Briefly, following minimization and a 200 ps equilibration, 20 ns simulations in water and in methanol were performed. The rotational autocorrelation times were calculated using a vector defined in the long dimension along the plane of the rings.

The Fab, Fc and the full-length antibody were fully solvated with SPC water molecules. Approximately 25,500 water molecules were used to solvate the Fab, 33,700 for Fc and 135,000 water molecules were used to solvate the full-length antibody. Chloride or sodium atoms were added to neutralize the overall charge of the system. Octahedral periodic boundary conditions were used in each of the simulations. The electrostatic interactions were calculated using PME with real space electrostatic cut off of 1.0 nm. The Lennard–Jones potential, describing the van der Walls interaction, was cut off at 1.0 nm. The Settle algorithm was used to constrain the bond lengths and angles of the water molecules; Lincs was used to constrain all other bond lengths, and the site algorithm, in Gromacs 4.0, was used to remove alkyl and amide hydrogen motions, allowing for the use of a time-step of 4 fs.

Throughout these simulations, the temperature was kept constant by coupling the system to a temperature bath of 300 K using the V-rescale algorithm. During equilibration, the pressure was kept constant by coupling the system to a pressure bath at 1.0 atm. Following equilibration, the simulations were kept at constant volume.

In all cases, following energy minimization, a 200 ps equilibration, at constant pressure, was performed to allow for the density of the system to converge. Following equilibration, a total of 24 trajectories were initiated and analyzed. Of the 24 trajectories, 19 were for the full-length antibody, four were for the Fab fragment and one was for the Fc fragment.

The full-length antibody trajectories were further broken down into two trajectories starting from a homology model based on 1HZH, two starting from homology models based on 1IGY, eight trajectories (4 1HZH based and 4 1IGY based) using structures taken from the ANM calculations (i.e., configurations along the modes calculated using ANM) and six trajectories in 150 mM sodium chloride.

The initial starting structures for the 150 mM salt simulation were selected from the collection of 1HZH and 1IGY trajectories. Three of the initial structures for the 150 mM salt calculations had no domain–domain interaction as defined as not having a salt bridge interaction of less than 5 Å. The other three initial structures had either a Fab1–Fab2 interaction, Fab1/Fc or Fab2/Fc interaction present as defined by at least one interdomain salt bridge having a distance of less than 5 Å.

All of full-length trastuzumab trajectories except for one 1HZH based trajectory, were simulated for approximately 80 ns. One of the 1HZH based trajectories was simulated for 500 ns. The 150 mM salt trajectories were simulated for 167 ns each. This resulted in a total simulation time of over 2 μs for the full-length antibody. The two dye conjugated Fab fragment simulations and one Fc simulation were each simulated for 400 ns. In two of the simulations, the Fab was conjugated with N-(1-pyrene)maleimide at LC and 205 and HC residue 121. The Fc was conjugated with N-(1-pyrene)maleimide at both HC 403 residues. Two additional Fab simulations without conjugated dyes were simulated for 200 ns and used to compare with the conjugated simulation. The complete list of the trajectories is given in supplementary Table 1.

The fluorescence anisotropy of the N-(1-pyrene)maleimide, attached to the Fab at LC 205 and HC residue 121 as well as to the Fc at HC 403, was calculated using methods described in Schroder et al., 2005. Briefly, the fluorescence anisotropy is given by

where μ is a vector representing both the normalized absorption and emission dipole moment vectors of the molecule. For the purposes of the calculation, the absorption and emission dipole moment vectors are estimated to be equivalent. For the case of N-(1-pyrene)-maleimide, the dipole moment vector used to calculate the fluorescence anisotropy is shown in supplemental Figure 3. All the anisotropy calculations were performed using the grotacf utility in Gromacs 4.0.

Solvent accessible surface area (SASA)

SASA was calculated using the g_sas code from gromacs

Mutual information

Inter-atomic mutual information (MI) for a trajectory from molecular dynamics simulation is defined as

where

is the displacement of the atomic position for the i-th atom relative to its average position over an ensemble; ρ(x) is the probability density of finding the i-th atom with a displacement of x; and ρ(x) is the probability density of finding the i-th atom with a displacement of x and the jth atom with a displacement of x. In our calculations we divided the displacements into discrete bins {Δr} and applied a discrete form of the equation (1):

The displacement bins for the atomic positions of the i-th atom are as follows: p(Δr) is the probability of finding the i-th atom with a displacement in the i-th bin Δr; p(Δr,Δr) is the probability of finding the i-th atom with a displacement in Δr and the j-th atom with a displacement in Δr. The summation is over all possible displacement bin combinations {Δr,Δr}. For each dimension, we tested the effects of number of bins and found that the patterns of correlations from MI do not change with bin numbers greater than 100. Hence, we chose 100 bins for each dimension.

Phi-psi entropy

For each non-terminal residue, the values for phi-psi angles were collected over all frames of the entire trajectory. For this phi-psi distribution, the phi-psi (φ-ψ) entropy is defined by the standard Shannon entropy:

Where p(φ,ψ) is the probability of the phi-psi angles falling in the { φ,ψ} bin. We used 3° interval to define bins in both phi and psi dimensions.

Inter-domain salt bridges

Trajectories with only the charged residues were extracted from full trajectories. The extracted trajectories were then imported into VMD, and the salt bridge analysis tools in VMD were used to dump the time courses for the salt bridges. We applied two cutoffs, 3.2 Å (VMD’s definition) and 5 Å to define the salt bridges. We then sorted the inter-domain salt bridges and calculated the various averages.

MSM construction

Choosing a lag time

Accurate construction of an MSM depends on choosing a set of conformational states in which the dynamics observed in simulations is sufficiently Markovian: i.e., the transition rate Tjk from state j to k (in some lag time τ) is independent of preceding transitions from states i to j (Tijk = TijTjk). Thus, the lag time τ, used to estimate transition rate is a key parameter. For example, a very short lag time may not sufficiently be able to account for the time needed for equilibration within each metastable state. As others have done previously,, we determined that τ = 1 ns is a suitable lag time by constructing a series of MSMs using a different lag times (1, 2 and 5 ns) and demonstrating that the rate of the slowest relaxation mode of the MSM is relatively insensitive to the choice of lag time (data not shown).

Solution of the master equation

The master equation describing the MSM dynamics is dp/dt = pK, where p is the state vector of populations, and K is an 8 × 8 matrix of rate coefficients. Equivalently, this can be expressed in terms of the transition matrix T = exp(τK), using the discrete propagation operation p(t + τ) = p(t)T, where Tij is the probability of transitioning from state i to state j in time τ. The solution of the master equation is p(t) = p(t = 0) B–1 exp(–Λt), where Λ is a matrix with the eigenvalues of K, λi, as diagonal entries, and B is the matrix of the (left) eigenvectors of K as columns. The properties of K are such that the largest eigenvalue is zero (corresponding to the stationary state, for which dp/dt = pK = 0), while the remaining eigenvalues λi are negative. The master equation kinetics can thus be described as a superposition of exponential relaxations for implied timescales τ*i = -λi−1. The sign structure (positive or negative entries) in each eigenvector describes the population flux occurring at each corresponding timescale. The second eigenvector (i.e., the largest non-zero λi) is the slowest kinetic relaxation in the system, and its corresponding eigenvector represents states exchanging population on this timescale. In practice, we calculate the implied timescales from the transition matrix T, using the relation

where µi are the eigenvalues of T.

Estimating the transition rate matrix

The entries of the rate matrix T are estimated from counts of transitions observed in the simulation. Finite sampling error in estimating transition rates is propagated as estimation error in the implied timescales, so it is this error analysis we are most interested in. Because the rate matrix must obey detailed balance (i.e., pi Tij = pj Tji, where pi are the equilibrium populations), inferring the transition matrix T from observed transition counts is a non-trivial problem, and several different methods exist, each with the purpose of inferring the posterior distribution of possible T given a finite collection of observed data., Here, we deal with this estimation problem in two ways. First, to ensure that our estimated rate matrix obeys detailed balance, we count all transitions i→j in both the forward and backward direction (this essentially assumes the simulation data are fully equilibrated, as molecular dynamics is microscopically reversible). Second, we estimate the posterior distribution of T nonparametrically, using the bootstrap method. Each bootstrap consists of drawing n trajectories, with replacement, from our finite set of n = 12 trajectories, and compiling the observed counts. Because trastuzumab is structurally symmetric (i.e., Fab1 = Fab2), Fab1–Fc and Fab2–Fc transitions were counted as the same. This bootstrap process was repeated 10,000 times to construct a posterior distribution of transition count matrices.

While the above bootstrap procedure addressed finite sampling error across trajectories, there is also sampling error due to the finite number of transitions observed. In particular, there are transitions (for example, to the three-contact state) that are not observed in the simulation data. To estimate entries of T for which no transitions are observed, we followed a pseudocount procedure in which, at each bootstrap iteration, instead of compiling the raw transition counts nij from each bootstrapped trajectory, the counts in each column j were drawn from a Dirichlet distribution

with parameters ui = nij + 1/N, where n = 8 is the number of states in the MSM. This expression is derived from the Bayesian posterior of a multinomial process given the observed transition counts, assuming a uniform prior.

Experimental

Protein expression and purification

Trastuzumab derivatives with engineered cysteine residues (THIOMABS) expressed and purified from Chinese Hamster Ovary (CHO) cells have been described previously. For E.coli expression of trastuzumab and trastuzumab “hingeless” variants, heavy and light chains were cloned into an expression vector and expressed as described., Proteins were purified by first capturing the antibody on protein A column (GE Healthcare, MAbSelect SuRe) and subsequently purifying by ion-exchange chromatography (GE Healthcare, HiTrap SP FF).

Site-specific conjugation of THIOMABs

Conjugation of the THIOMABs at engineered cysteines has been described previously. Briefly, THIOMABs were reduced with 20-fold molar excess of DTT at pH 7.0, 25°C over 14–16 h. Reducing agent was removed by desalting (GE Healthcare, PD-10). Native disulfide bridges were reformed by mild re-oxidation using 10-fold molar excess dhAA (Sigma Aldrich) over 3 h at 25°C. Conjugation reaction was initiated by addition of 3-fold molar excess (complete labeling) of N-(1-Pyrene)maleimide (Invitrogen). Partial labeling was achieved using equal molar amount of N-(1-Pyrene)maleimide (Invitrogen). Conjugation was performed at 25°C for 30 min. Any unreacted cysteine residues were alkylated using 10-fold molar excess of sodium iodoacetate (Sigma Aldrich). Excess dye was removed by desalting (GE Healthcare, PD-10). Conjugation efficiency was determined by measuring protein and dye absorbance (at 280 nm and 342 nm, respectively). Samples were concentrated and stored in 50 mM sodium acetate buffer pH5.5 at 5°C. Conjugation specificity was confirmed by mass spectrometry.

Size-exclusion chromatography

Unconjugated trastuzumab variants were analyzed using TSKgel G3000SWXL column (10 μm, 7.8 mm × 30 cm, TOSOH Biosciences) on Agilent 1100 HPLC system (operated using Chemstation software package). Mobile phase used: 0.25 M potassium phosphate buffer pH 6.9, 0.5 M potassium chloride. The flow rate was 0.5ml/min. These conditions were used for analysis of E.coli-derived hinged and “hingeless” THIOMAB LC205 and for monitoring the formation of intermolecular disulfide bonds during the re-oxidation step of THIOMAB conjugation.

Pyrene conjugation causes significant interaction of the antibody with the above gel filtration resin. For analysis of conjugated molecules, we chose to use non-silica resin: and Superdex 200 10/30 GL (GE Healthcare) with PBS as a mobile phase at flow rate 0.5 ml/min.

Hydrophobic-interaction chromatography (HIC)

Unconjugated THIOMAB and THIOMAB conjugates with single or two pyrene probes were separated on TSKgel Phenyl NPR-5PW column (Tosoh Biosciences, 10 μm, 7.5 mm × 7.5 cm) using linear gradient from 0–35%B over 200 min (buffer A: 1M ammonium sulfate, 50 mM sodium phosphate pH 7.0, buffer B: 50 mM sodium phosphate, 25% isopropanol) at flow rate 0.8 ml/min. Up to 20 mg of protein sample was loaded and purified on ÄKTA FPLC system (GE Healthcare).

Fluorescence anisotropy measurements

Sample preparation

Pyrene-labeled trastuzumab THIOMABs were diluted from the stock solution (100 μM) into appropriate buffer to final 1–3 μM in 1 ml quartz cuvette (Starna Cells) for fluorescence measurements. Trastuzumab molecules with reduced disulfide bonds were prepared by treating the antibody solution with 20× molar excess DTT for an appropriate amount of time (20 min or 12 h) at pH 7.0 at room temperature and alkylating the exposed disulfides with 100× molar excess of sodium iodoacetate. Samples were desalted (GE Healthcare, PD-10) and stored in 50 mM sodium acetate buffer pH 5.5 at 5°C.

Data acquisition

Time-resolved anisotropy was measured using time-correlated single-phonon counting technique (TCSPC) using FluoroMax4 (HORIBA Jobin Yvon) steady-state spectrofluorometer equipped with FluoroLog TCSPC lifetime module. Data was collected using Data Station software supplied by the manufacturer. A 340 nm NanoLED was used as an excitation source (with optical pulse duration of 1.4 ns), and decays were collected at 415 nm at 10 nm slit width. Lifetime decays were collected with excitation and detection polarizers at magic angle conditions. For anisotropy decays, two time decays were measured, parallel (I) and perpendicular (I) (subscripts denote orientations of excitation and emission polarizers, respectively; V, vertical; H, horizontal). The correction factor (G factor) was determined by measuring I and I. Instrument response function (IRF) was recorded using a dilute solution of colloidal silica (Sigma Aldrich). Decay data was analyzed using DAS6 software package (HORIBA Jobin Yvon). Briefly, the anisotropy decay law r(t):

was determined using impulse reconvolution analysis by first fitting the denominator S(t) (total fluorescence decay) to the multiexponential decay using the method of nonlinear least squares. Using the obtained total fluorescence decay, the product r(t)S(t) was convoluted with the IRF to find the best fit to the experimental difference decay D(t). A single- or two-exponential anisotropy decay was used for impulse reconvolution analysis:

Average anisotropy decay time

was used as a measure of overall flexibility of the full-length antibody or the antibody fragment with a fluorescence probe at a given conjugation position.

Conclusions

We used experimental and computational methods to characterize the solution dynamics of the IgG1 antibody trastuzumab and find that the full-length antibody exists as a heterogeneous population of meta-stable conformations with non-specific domain–domain interactions between fragments. We measured the hydrodynamic parameters of fragments and full-length antibody and showed that the same parameters estimated from all-atom molecular dynamics simulations were quantitatively similar, thus validating the computational models. We further analyzed the atomic details of our simulation data by looking at correlated motion between residue pairs and phi/psi conformational entropy. Dynamics in our models are consistent with previous data such as high structural flexibility in the hinge region and heavy chain CDR-3. Correlated motion suggests that Ig domains are highly correlated along the direction of their β strands (e.g., heavy chain variable domain to light chain variable domain) and less correlated in the “end-to-end” direction (e.g., variable domain to CH1). This result indicates that domains that are functionally coupled (e.g., variable domains that bind to a single antigen) have coupled dynamics, and the dual functions of the antibody (antigen binding and effector binding) are dynamically independent. Importantly, we show significant and non-specific inter-fragment domain–domain interactions. We incorporate the molecular dynamics data (> 1 μs, an order of magnitude more than what was previously published) into a Markov Model and show an equilibrium between multiple meta-stable conformations of full-length trastuzumab with domain–domain interactions between fragments. The biological significance of these domain–domain interactions remains to be determined.