Kinetics and Energetics of Biomolecular Folding and Binding

The ability of biomolecules to fold and to bind to other molecules is fundamental to virtually every living process. Advanced experimental techniques can now reveal how single biomolecules fold or bind against mechanical force, with the force serving as both the regulator and the probe of folding and binding transitions. Here, we present analytical expressions suitable for fitting the major experimental outputs from such experiments to enable their analysis and interpretation. The fit yields the key determinants of the folding and binding processes: the intrinsic on-rate and the location and height of the activation barrier.Dynamic processes in living cells are regulated through conformational changes in biomolecules—their folding into a particular shape or binding to selected partners. The ability of biomolecules to fold and to bind enables them to act as switches, assembly factors, pumps, or force- and displacement-generating motors (1). Folding and binding transitions are often hindered by a free energy barrier. Overcoming the barrier requires energy-demanding rearrangements such as displacing water from the sites of native contacts and breaking nonnative electrostatic contacts, as well as loss of configurational entropy. Once the barrier is crossed, the folded and bound states are stabilized by short-range interactions: hydrogen bonds, favorable hydrophobic effects, and electrostatic and van der Waals attractions (2).Mechanistic information about folding and binding processes is detailed in the folding and binding trajectories of individual molecules: observing an ensemble of molecules may obscure the inherent heterogeneity of these processes. Single-molecule trajectories can be induced, and monitored, by applying force to unfold/unbind a molecule and then relaxing the force until folding or binding is observed (3–5) (Fig. 1). Varying the force relaxation rate shifts the range of forces at which folding or binding occurs, thus broadening the explorable spectrum of molecular responses to force and revealing conformational changes that are otherwise too fast to detect. The measured force-dependent kinetics elucidates the role of force in physiological processes (6) and provides ways to control the timescales, and even the fate, of these processes. The force-dependent data also provides a route to understanding folding and binding in the absence of force—by extrapolating the data to zero force via a fit to a theory.Open in a separate windowFigure 1Schematic of the output from a force-relaxation experiment. The applied force is continuously relaxed from the initial value F₀ until the biomolecule folds or binds, as signified by a sharp increase in the measured force. From multiple repeats of this experiment, distributions of the folding or binding forces are collected (inset). Fitting the force distributions with the derived analytical expression yields the key parameters that determine the kinetics and energetics of folding or binding.In this letter, we derive an analytical expression for the distribution of transition forces, the major output of force-relaxation experiments that probe folding and binding processes. The expression extracts the key determinants of these processes: the on-rate and activation barrier in the absence of force. The theory is first developed in the context of biomolecular folding, and is then extended to cover the binding of a ligand tethered to a receptor. In contrast to unfolding and unbinding, the reverse processes of folding and binding require a theory that accounts for the compliance of the unfolded state, as well as the effect of the tether, to recover the true kinetic parameters of the biomolecule of interest.In a force-relaxation experiment, an unfolded biomolecule or unbound ligand-receptor complex is subject to a stretching force, which is decreased from the initial value F₀ as the pulling device approaches the sample at speed V until a folding or binding transition is observed (Fig. 1) (3–5). Define S(t) as the probability that the molecule has not yet escaped from the unfolded (implied: or unbound) state at time t. When escape is limited by one dominant barrier, S(t) follows the first-order rate equationS˙(t)≡dS(t)dt=−k←(F(t))S(t),where k_←(F(t)) is the on-rate at force F at time t. Because, prior to the transition, the applied force decreases monotonically with time, the distribution of transition forces, p(F), is related to S(t) through p(F)dF=S˙(t)dt, yieldingp(F)=−k←(F)F˙(F)e−∫F0Fk←(F′)F˙(F′)dF′.(1)Here F˙(F)≡dF(t)/dt<0 is the force relaxation rate. The proper normalization of p(F) is readily confirmed by integrating Eq. 1 from the initial force F₀ to negative infinity, the latter accounting for transitions that do not occur by the end of the experiment. Note that the expression for the distribution of folding/binding forces in Eq. 1 differs from its analog for the unfolding process (7) by the limits of integration and a negative sign, reflecting the property of a relaxation experiment to decrease the survival probability S(t) by decreasing the force. Converting the formal expression in Eq. 1 into a form suitable for fitting experimental data requires establishing functional forms for k_←(F) and F˙(F) and analytically solving the integral. These steps are accomplished below.The on-rate k_←(F) is computed by treating the conformational dynamics of the molecule as a random walk on the combined free energy profile G(x,t) = G₀(x) + G_pull(x,t) along the molecular extension x. Here G₀(x) is the intrinsic molecular potential and G_pull(x,t) is the potential of the pulling device. When G(x,t) features a high barrier on the scale of k_BT (k_B is the Boltzmann constant and T the temperature), the dynamics can be treated as diffusive. The unfolded region of the intrinsic potential for a folding process, unlike that for a barrierless process (8), can be captured by the functionG0(x)=ΔG‡ν1−ν(xx‡)11−ν−ΔG‡ν(xx‡),which has a sharp (if ν = 1/2, Fig. 2, inset) or smooth (if ν = 2/3) barrier of height ΔG^‡ and location x^‡. The potential of a pulling device of stiffness κ_S is G_pull(x,t) = κ_S/2(X₀ – Vt – x)² with an initial minimum at X₀ (corresponding to F₀). Applying Kramers formalism (9) to the combined potential G(x,t), we establish the analytical form of the on-rate at force F(t),k←(F)=k0(1+κSκU(F))1ν−12(1+νFx‡ΔG‡)1ν−1×eβΔG‡[1−(1+κSκU(F))2ν1−ν−1(1+νFx‡ΔG‡)1ν],where k₀ is the intrinsic on-rate, β ≡ (k_BT)⁻¹, andκU(F)=ν(1−ν)2ΔG‡x‡2(1+νFx‡ΔG‡)2−1νis the stiffness of the unfolded biomolecule under force F (see the Supporting Material for details on all derivations). The full nonlinear form of G_pull(x,t) was necessary in the derivation because, in contrast to the typically stiff folded state, the unfolded state may be soft (to be exact, 1/2κ_S x^‡2(F) << k_BT may not be satisfied) and thus easily deformed by the pulling device. Because of this deformation, the folding transition faces an extra contribution (regulated by the ratio κ_S/κ_U(F)) to the barrier height, typically negligible for unfolding, that decreases the on-rate in addition to the applied force F.Open in a separate windowFigure 2Contributions to the free energy profile for folding (inset) and binding (main figure). The derived expression (Eq. 2) extracts the on-rate and the location and height of the activation barrier to folding. When applied to binding data, the expression extracts the parameters of the ligand-tether-receptor (LTR) potential G˜0 (x); the proposed algorithm (Eqs. 3 and 4) removes the contribution of the tether potential G_teth(x) to recover the parameters of the intrinsic ligand-receptor (LR) potential G₀(x).The last piece required for Eq. 1, the loading rate F˙(F), is computed as the time derivative of the force F(t) on the unfolded molecule at its most probable extension at time t:F˙(F)=−κSV1+κS/κU(F).Finally, we realize that the integral in Eq. 1 can be solved analytically exactly, both for ν = 1/2 and ν = 2/3, resulting in the analytical expression for the distribution of folding forces:p(F)=k←(F)|F˙(F)|e−k←(F)β|F˙(F)|x‡(1+κSκU(F))νν−1(1+νFx‡ΔG‡)1−1ν.(2)Equation 2 can be readily applied to (normalized) histograms from force-relaxation experiments to extract the parameters of the intrinsic kinetics and energetics of folding. Being exact for ν = 1/2 and ν = 2/3, Eq. 2 is also an accurate approximation for any ν in the interval 1/2 < ν < 2/3 as long as κ_S ≲ κ_U (F) (see Fig. S1 in the Supporting Material). For simplicity, in Eq. 2 we have omitted the term containing F₀ as negligible if F₀ is large enough to prevent folding events.The solution in Eq. 2 reveals properties of the distribution of folding forces that distinguish it from its unfolding counterpart (7):

• 1.The distribution has a positive skew (Fig. 3), as intuitively expected: the rare folding events occur at high forces when the barrier is still high.Open in a separate windowFigure 3Force histograms from folding (left) and binding (right) simulations at several values of the force-relaxation speed (in nanometers per second, indicated at each histogram). Fitting the histograms with the analytical expression in Eq. 2 (lines) recovers the on-rate and activation barrier for folding or binding (2.Increasing the relaxation speed shifts the distribution to lower forces (Fig. 3): faster force relaxation leaves less time for thermal fluctuations to push the system over a high barrier, causing transitions to occur later (i.e., at lower forces), when the barrier is lower.
• 3.The stiffness κ_S and speed V enter Eq. 2 separately, providing independent routes to control the range of folding forces and thus enhance the robustness of a fit.

The application of the above framework to binding experiments on a ligand and receptor connected by a tether (3) involves an additional step—decoupling the effect of the tether—to reconstruct the parameters of ligand-receptor binding. Indeed, the parameters extracted from a fit of experimental histograms to Eq. 2 characterize the ligand-tether-receptor (LTR) potential (k˜0, x˜‡, ΔG˜‡, ν) (Fig. 2). The parameters of the natural ligand-receptor (LR) potential (k₀, x^‡, ΔG^‡) can be recovered using three characteristics of the tether: contour length L; persistence length p; and extension Δℓ of the tether along the direction of the force in the LTR transition state. The values of L and p can be determined from the force-extension curve of the tether (10); these define the tether potential G_teth(x) (Fig. 2). The value of Δℓ can be found from an unbinding experiment (7) on LTR and the geometry of the tether attachment points (see Fig. S3). Approximating the region of the LR potential between the transition and unbound states as harmonic, with no assumptions about the shape of the potential beyond x^‡, the ligand-receptor barrier parameters are thenx‡=α−1α−2x˜‡,ΔG‡=(α−1)22(α−2)x˜‡Fteth(Δℓ+x˜‡),(3)and the intrinsic unimolecular association rate isk0≈k˜0(βΔG‡)32(βΔG˜‡)1ν−12(x˜‡x‡)2eβ(ΔG˜‡−ΔG‡).(4)Here, the force value Fteth(Δℓ+x˜‡) is extracted from the force-extension curve of the tether at extension Δℓ+x˜‡ andα=2(ΔG˜‡−Gteth(Δℓ)+Gteth(Δℓ+x˜‡))x˜‡Fteth(Δℓ+x˜‡),where G_teth(x) is the wormlike-chain potential (see Eq. S13 in the Supporting Material). Equations 3–4 confirm that a tether decreases the height and width of the barrier (see Fig. 2), thus increasing the on-rate.In Fig. 3, the developed analytical framework is applied to folding and binding force histograms from Brownian dynamics simulations at parameters similar to those in the analogous experimental and computational studies (3,5,11) (for details on simulations and fitting procedure, see the Supporting Material). For the stringency of the test, the simulations account for the wormlike-chain nature of the molecular unfolded and LTR unbound states that is not explicitly accounted for in the theory. With optimized binning (12) of the histograms and a least-squares fit, Eqs. 2–4 recover the on-rate, the location and the height of the activation barrier, and the value of ν that best captures how the kinetics scale with force (

1.Multiple relaxation speeds,

2.Folding/binding events at low forces, and

3.A large number of events at each speed.

Table 1

On-rate and the location and height of the activation barrier from the fit of simulated data to the theory in Eq. 2

The Physical Foundation of Vasoocclusion in Sickle Cell Disease

The pathology of sickle cell disease arises from the occlusion of small blood vessels because of polymerization of the sickle hemoglobin within the red cells. We present measurements using a microfluidic method we have developed to determine the pressure required to eject individual red cells from a capillary-sized channel after the cell has sickled. We find that the maximum pressure is only ∼100 Pa, much smaller than typically found in the microcirculation. This explains why experiments using animal models have not observed occlusion beginning in capillaries. The magnitude of the pressure and its dependence on intracellular concentration are both well described as consequences of sickle hemoglobin polymerization acting as a Brownian ratchet. Given the recently determined stiffness of sickle hemoglobin gels, the observed obstruction seen in sickle cell disease as mediated by adherent cells can now be rationalized, and surprisingly suggests a window of maximum vulnerability during circulation of sickle cells.Human capillaries are narrower than the erythrocytes they convey. In sickle cell disease, red cells can become rigid in those capillaries, because the hemoglobin inside the red cell will aggregate into stiff polymers. This happens once the molecules deliver their oxygen, and led to the long-held view that capillary occlusion was central to the pathophysiology of the disease (1,2). This was challenged when microscopic study of animal model tissues perfused with sickle blood revealed blockages that began further downstream, in the somewhat larger venules (3–5), at the site of adherent red or white cells which diminished the vessel lumen without fully obstructing the flow. Yet no rationale has been presented for the failure of the prior assumption of capillary blockage. Microfluidic methods (6) are ideally suited to discover why cells don’t get stuck in the capillaries, yet occlude subsequent vessels, and we have constructed a system to address this question. Our measurements show that the pressure differences across capillaries in vivo can easily dislodge a cell sickled within a capillary, giving an experimental answer to the question of why sickled cells don’t stick in capillaries. It turns out that the pressure a cell can withstand is quantitatively explained by the Brownian ratchet behavior of sickle hemoglobin polymerization.We constructed single-cell channels in transparent polydimethylsiloxane, with a cross section (1.5 μm × 4 μm) that is smaller than the resting diameter of red cells (Fig. 1). These channels are much narrower than those that have been employed in other recent studies of the sickling process (7,8), and they resemble human capillaries in permitting only one cell at a time to pass through them. We used a laser photolysis method to create ligand free (deoxygenated) cells, and this requires that the hemoglobin bind CO, which can then be readily removed by strong illumination, in contrast to bound O₂ which is released with far lower efficiency than CO. The microfluidic chips were enclosed in a gas-tight chamber flushed with CO to avoid introduction of oxygen and keep the cells fully ligated before photolysis. The profiles of the channels were confirmed by microscopic observation. To confirm that liquid did not pass around the cells when they were trapped in the channels, fluorescent beads were introduced into some cell solutions. The beads did not pass the cells, nor did they approach the cell when it was occluded, verifying that no significant flow occurred around the cell when it was stuck.Open in a separate windowFigure 1An erythrocyte enters a channel (moving left to right) and is positioned in the center, where it will be photolyzed. The channel cross section is 1.5 μm × 4 μm, smaller than a resting red cell diameter.Optical measurements were carried out on a microspectrophotometer constructed on an optical table. The system employed ×32 LWD objectives (Leitz, Wetzlar, Germany), which were autofocused during collection of absorption spectra to minimize aberrations. Spectra were obtained using a series 300 camera (Photometrics, Tucson, AZ); video imaging was done with a high-speed camera (Photron, San Diego, CA). Photolysis was provided by a 2020 Argon Ion laser (Spectra Physics, Houston, TX). Sickle cells were obtained from patients at the Marian Anderson Sickle Cell Center at St. Christopher''s Hospital for Children, Philadelphia, PA by phlebotomy into EDTA-containing tubes. The blood was centrifuged at 5°C at 1200g for 4 min, and then the pellet was washed 4× with 15 volumes of buffer (120 mM NaCl, 2 mM KCl, 10 mM dibasic Na Phosphate, 7 mM monobasic Na Phosphate, 3.4 mM Na Bicarbonate, and 6 mM Dextrose) by repeated suspension and centrifugation at 30g for 4 min. This minimizes fibrinogen and platelets in the final suspension, to insure that these studies are controlled by the mechanical properties of the cells themselves.Our experiment began by parking a cell in the center of a channel (Fig. 1). The cell, its hemoglobin, and the microchannel environment all were saturated with CO. Because the thickness of the channel is known, we were able to determine the hemoglobin concentration inside the cell from its absorption spectrum (Fig. 2 A). Steady-state laser illumination then removed the CO, allowing the hemoglobin to polymerize, in which condition it remained while the laser was kept on. Removal of CO was confirmed by observing the spectral difference between COHb and deoxyHb. Photolysis of COHb generates negligible heating (9–11). During illumination, hydrostatic pressure was applied until the cell broke free.Open in a separate windowFigure 2(A) Absorption of the cell (points), fit to a standard spectrum (9). (B) Pressure to dislodge a cell sickled in the microchannel, as a function of intracellular concentration. Note that typical intracellular concentrations are ∼32 g/dL. (Line) Brownian-ratchet theory described in the text. The coefficient of friction (0.036) is within the observed range, and is the only parameter varied.The magnitude of the dislodging pressure, measured by simple height difference between input and output cell reservoirs, is shown in Fig. 2 B. The pressure needed to dislodge the cell increased with increasing intracellular Hb concentration, implying that an increased mass of polymerized hemoglobin is more difficult to dislodge. A clear concentration threshold for capture is apparent. While there is a well-known solubility below which polymers cannot form (18.5 g/dL for the 22°C of this experiment (12)), the threshold here is significantly higher.Central to explaining these observations is a Brownian ratchet mechanism (13) which derives from the metastable nature of this polymerization process. Unless disrupted, as by centrifugation, polymerization in sickle hemoglobin terminates before the thermodynamic limit of monomer solubility is reached (14,15). This arises from the fact that polymers only grow at their ends, which are easily occluded in the dense mass of polymers that form. This end obstruction leaves the system in a metastable state and fluctuations accordingly provide polymers with space into which they can incrementally grow. This Brownian ratchet has been shown to lead to dramatic fiber buckling when individual fibers are isolated in sickle cells (16). The force can be simply expressed as f = (kT/δ) ln S(c), where k is Boltzmann’s constant, T the absolute temperature, δ the net spatial elongation from addition of a single monomer, and S is the supersaturation of the solution when the metastable limit is reached, at monomer concentration c. In this calculation, c is taken as the terminal concentration, computed from our empirical finding (15) that in this metastable system the amount of polymerized hemoglobin Δ is Δ(∞) = 2/3 (c_o-c_s), rather than the expected thermodynamic limit c_o-c_s, where c_o is the initial concentration and c_s is the solubility.For determining the net force, the total number of fibers must be known, and can be calculated based on the double nucleation mechanism (17) which has been quantitatively successful in describing polymerization. The concentration of polymers [p(t)] initially grows exponentially, described by[p(t)]=(AB2J)exp(Bt),where A and B are parameters related to nucleation, and J is the polymer elongation rate, as described in Ferrone et al. (17). Because A and B are both extremely concentration-dependent (9), they will drop dramatically once monomers begin to add to polymers in any significant numbers, and thereby diminish the remaining monomer pool. Thanks to the extreme concentration dependence of the reaction, this rapidly shuts off further polymerization. This happens at approximately the 10th time (the time when the reaction has reached 1/10 of its maximum). Thus, the [p(t_1/10)] ≈ [p(∞)]. Moreover, at one-tenth of the reaction,Δ(t1/10)=12Aexp(Bt1/10)=Δ(∞)10,and thus[p(∞)]=(BJ)(Δ(∞)10)=(BJ)((co−cs)15).For computing the number of fibers, the volume of the cell was taken as 90 μm³. This calculation shows, as expected, that the number of polymers in the cell is highly concentration-dependent, and very few fibers are produced at concentrations just above solubility, but the number grows sharply as concentration rises. This is the main contribution to the threshold in holding force shown by the data.With the force per fiber, and the total number of fibers, the net force against the wall is known. With a coefficient of friction, this reveals the force that a trapped cell can withstand. If the force is divided by the cross-sectional area across which the force is applied, we get a prediction of the dislodging pressure, which can be compared to the data. For a quantitative comparison with the results, two further corrections, of order unity, were applied. Because only normal force will contribute to friction, the calculated force was determined by integrating cos θ. This integration is not over all angles (π) because of the possibility that large incidence angles of the fibers against the wall will lead to fiber runaway (18). Therefore, the integration described is taken to the runaway threshold, here ∼1 rad. Finally, it is necessary to assign a coefficient of friction. Known values span the range of 0.03–0.06 (19). We therefore selected a value within the range, 0.036, as the best match for the data. The predicted pressures match the measurements well, as the line in Fig. 2 B shows.Because the flow resistance is comparable for red cells traversing glass channels and endothelial-lined capillaries (20), we conclude that in vivo the pressures a sickled cell inside a capillary can withstand are no more than hundreds of Pa. This is significantly smaller than typical arteriovenous pressure differentials that have been measured, which range from 0.7 kPa (in hamster skin (21)) to 7.9 kPa (in rat mesentery (22)).Our measurements coupled with recent determination of the stiffness of sickle hemoglobin gels (23) provide the missing physical basis for the processes of vasoocclusion seen in ex vivo tissue and animal models of sickle cell disease, arguing that these observations indeed represent fundamental behavior of sickle cell disease. We now understand this behavior in terms of three possible outcomes, all intimately connected with kinetics:

• 1.Certain escape: A cell that does not polymerize until after passing the obstruction can reach the lungs where it reoxygenates and resets its polymerization clock.
• 2.Possible escape: A cell that polymerizes within the capillary will assume an elongated sausage shape. The forces that it can exert against the wall cannot hold it there, and it will emerge into the postcapillary venule. There it has some chance of passing a subsequent obstruction, though it might also obstruct flow were it to rotate before reaching the adherent cell, so as to present its long dimension to the reduced space it must traverse.
• 3.Certain occlusion: A cell that does not polymerize in the capillary reassumes a larger diameter as soon as it escapes. If the cell then polymerizes before it encounters a cell attached to the venule wall, this rigidified cell will not be able to squeeze past the adherent cell, because that kind of deformation takes MPa (23). This would precipitate the type of blockage that is observed. This suggests that there is a window of greatest vulnerability, toward which therapies might be addressed.

     

Quantitative Analysis of Three-Dimensional Fluorescence Localization Microscopy Data

Identifying the three-dimensional molecular organization of subcellular organelles in intact cells has been challenging to date. Here we present an analysis approach for three-dimensional localization microscopy that can not only identify subcellular objects below the diffraction limit but also quantify their shape and volume. This approach is particularly useful to map the topography of the plasma membrane and measure protein distribution within an undulating membrane.Single molecule localization microscopy (SMLM) (1–3) is a superresolution fluorescence microscopy technique that produces coordinate data for single molecule localizations with a precision of tens of nanometers in live and fixed cells. These methods have mainly been performed with total internal reflectance fluorescence microscopy and therefore have generated two-dimensional molecular coordinates. Such two-dimensional data sets have revealed nanosized clusters of membrane proteins at the cell surface (4–7). This was achieved with analysis routines based on pair-correlation analysis (8), Ripley’s K function (9), and related techniques. While three-dimensional localization microscopy techniques such as biplane imaging (10), astigmatic spot analysis (11), and depth-encoding point-spread functions (12) have now been developed, quantitative analysis approaches of three-dimensional coordinate patterns have not.Here, we describe an approach based on Getis and Franklin''s local point pattern analysis to quantitatively analyze three-dimensional subcellular structures and map plasma membrane topography. The latter can also be used to account for topography-induced clustering of membrane proteins in an undulating membrane. To illustrate the approach, we generated three-dimensional SMLM data of the membrane dye DiI and the protein Linker for Activation of T cells (LAT) fused to the photoswitchable fluorescent protein mEos2 in T cells. It has been previously shown that LAT resides within the plasma membrane as well as membrane-proximal vesicles (5,13). The data were acquired using the biplane SMLM technique and highly inclined and laminated optical sheet illumination (14). Three-dimensional molecular coordinates were calculated by fitting a three-dimensional theoretical point-spread-function to the acquired data.As previously described for two-dimensional SMLM data analysis (5), Ripley’s K-function is calculated according to Eq. 1 where V is the analyzed volume, n is the total number of points, and r is the radius of a sphere (a circle for the two-dimensional case) centered on each point. The value K(r) is thus a measure of how many points are encircled within a sphere of radius r:K(r)=V∑i=1n∑j=1n(δij/n2);δij=1ifd(pointi,pointj)<r,0else.(1)For completely spatially random (CSR) data, K(r) scales with the volume of the sphere. We therefore linearize the K-function such that it scales with radius (the L-function) using:L(r)=(3K(r)4π)1/3.(2)The value of L(r)−r is then zero for the CSR case. Values of L(r)−r above zero indicate clustering at the length scale, r.Next we used the related Getis and Franklin''s local point pattern analysis to generate a clustering value (L(r) at r = 50 nm; L(50)) for each point, j, based on the local three-dimensional molecular density. This was calculated using:Lj(50)=((3V4π)∑i=1n(δijn))1/3;δij=1ifd(pointi,pointj)<50,0else.(3)These values can then be interpolated such that every voxel in a volume is assigned a cluster value based on the number of encircled points, relative to the expected CSR case. This allows construction of isosurfaces where all points on the surface have an identical L(50) value. A high threshold imparts a strict criterion for cluster detection compared to a lower one, and this allows users to, for example, determine the efficiency of sequestration into clusters by quantifying the cluster number and size as a function of the threshold.To illustrate the identification of subcellular structures, Lat-mEos2 was imaged by three-dimensional SMLM in activated T cells at the immunological synapse (Fig. 1 A). Three-dimensional projections of isosurfaces (for L(50) = 200) clearly identified intracellular LAT vesicles at varying depths within the synapse (Fig. 1, B and C). Cluster statistics were extracted from this data set to quantify the distribution of clusters in the z direction as well as the volume and sphericity of the LAT objects themselves (Fig. 1, D–F).Open in a separate windowFigure 1Identification of subcellular objects in three dimensions by isosurface rendering of molecular distribution. (A) Schematic of a T cell synapse formed against an activating coverslip where subsynaptic LAT vesicles (red dots) can be imaged with three-dimensional SMLM. (B and C) Isosurfaces, shown in x,z view (B) and as projection (C), identify T cell vesicles as LAT objects with L(50) > 200 (Eq. 3). (D–F) Distribution of LAT objects in z direction (D), volume (E), and sphericity (F) of LAT objects in T cells.Membrane undulations can cause clustering artifacts when the distribution of membrane proteins is recorded as a two-dimensional projection (15) (Fig. 2 A), as is the case in two-dimensional SMLM under total internal reflectance fluorescence illumination. To illustrate a solution to this problem, we obtained three-dimensional SMLM data sets of the membrane dye DiI (16) in resting T cells adhered onto nonactivating coverslips. With appropriately short labeling times to prevent dye internalization, it can be assumed that all DiI molecules reside in the plasma membrane. In this case, as is the case for plasma membrane proteins, neither two-dimensional nor three-dimensional analysis is appropriate, as it is a priori known the points must be derived from a two-dimensional membrane folded in three-dimensional space. To correct for membrane undulations, the plasma membrane topography must first be mapped so that molecular coordinates of membrane molecules can be appropriately corrected in two-dimensional projections. The position of the plasma membrane in three dimensions, i.e., the membrane topography, was determined by averaging the z position of all DiI molecules within a 100-nm radius in x-y at each point. The averaged z-position of DiI molecules was then displayed as a map, which exhibits a smooth, undulating profile (Fig. 2 B). The selection of this radius determines the accuracy of the assigned z position but also causes smoothing of the membrane profile.Open in a separate windowFigure 2Mapping of membrane topography and correction of molecular distributions in undulating membranes. (A) Two-dimensional projections can cause cluster artifacts, for example in membrane ruffles. Molecules (red rectangles) in the upper image are equally spaced along the membrane but appear as clusters in two-dimensional projections in areas with high gradient. (B) Three-dimensional membrane topography of a 2 × 2 μm plasma membrane area of a resting T cell obtained from averaged z positions of DiI molecules. Note that membrane undulation is ∼100 nm. (C) Map of membrane gradient, corresponding to the topography map shown in panel B, with an area of high gradient highlighted (dashed red box). (D) Correction of the circle radii in the Getis and Franklin cluster map calculations to account for projection artifacts. (E and F) Cluster map of data shown in panel C before (E) and after (F) correction for membrane gradient. Boxes in panels C, E, and F highlight the regions with high membrane gradient.Next, the gradient at the position of each DiI molecule was determined and interpolated into a gradient map (Fig. 2 C). Here, blue represents horizontal, i.e., flat membrane areas, whereas red regions indicate areas of high gradient. The information from the gradient map was then used to ensure that the two-dimensional circles in the Getis and Franklin cluster map calculations each correspond to an identical area of membrane, hence accounting for two-dimensional projection artifacts. To do this, the size of the circle (r) used to calculate the L value for each molecule was modified using Eq. 4, where c is calculated for the surface, S, using Eq. 5:r(corr)=r(uncorr)(1+c2)1/4,(4)c=((∂S∂x)2+(∂S∂y)2)1/2.(5)This operation is shown schematically in Fig. 2. The comparison of Getis and Franklin cluster maps before (Fig. 2 E) and after (Fig. 2 F) correction for the gradient shows that cluster values for DiI molecules were substantially reduced by up to 5–10% at sites where the plasma membrane had a high gradient (area highlighted in red box), and where the two-dimensional projection of three-dimensional structures caused an overestimation of clustering.In conclusion, we demonstrated that three-dimensional superresolution localization microscopy data can be used to identify and quantify subcellular structures. The approach has the distinct advantage that subcellular structures are solely identified by the distribution of the fluorescent marker so that no a priori knowledge of the structure is necessary. How precisely the subcellular structures are identified only depends on how efficiently the fluorescent maker is recruited to the structure, and hence does not depend on the resolution limits of optical microscopy. We applied the methods to two very different structures in T cells: small intracellular vesicles and the undulating plasma membrane. Importantly, the topography of plasma membrane can also be used to correct clustering artifacts in two-dimensional projections, which may be useful for distribution analysis within membranes.     

Indolicidin Binding Induces Thinning of a Lipid Bilayer

We use all-atom molecular dynamics simulations on a massive scale to compute the standard binding free energy of the 13-residue antimicrobial peptide indolicidin to a lipid bilayer. The analysis of statistical convergence reveals systematic sampling errors that correlate with reorganization of the bilayer on the microsecond timescale and persist throughout a total of 1.4 ms of sampling. Consistent with experimental observations, indolicidin induces membrane thinning, although the simulations significantly overestimate the lipophilicity of the peptide.Antimicrobial peptides are a component of the innate immune system of eukaryotes (1). As such, they must interact with pathogenic membranes, either during translocation or by disrupting their structural integrity (2). Here we examine the binding of the 13-residue cationic antimicrobial peptide indolicidin (3) (ILPWKWPWWPWRR-NH₂) to a lipid membrane as a first step towards elucidating its mechanism of action.Molecular solutes interact with lipid membranes in many cellular processes (4). Computational approaches such as molecular dynamics simulations have been widely used to characterize these interactions (5). However, molecular dynamics simulations can require unfeasibly long times to reach equilibrium (6). Therefore, it is common to compute equilibrium properties of solute insertion into lipid bilayers using umbrella sampling (7) simulations in which the solute is restrained along the bilayer normal using harmonic restraining potentials, or umbrellas, centered at zi0 values distributed between bulk water and the bilayer center.It is often assumed that equilibrium properties rapidly attain convergence in umbrella sampling simulations; accordingly, convergence measures are rarely published (8). However, we have recently shown that umbrella sampling simulations require up to 100 ns per umbrella (3 μs in total) to eliminate systematic sampling errors in the standard free energy of binding, ΔGbind0, of an arginine side-chain analog from bulk water to a lipid bilayer (8). The fact that umbrella sampling has been used to investigate the bilayer insertion of substantially larger solutes (9) motivates a systematic evaluation of statistical sampling convergence of ΔGbind0 for indolicidin in a lipid bilayer.To estimate ΔGbind0 of indolicidin to a lipid bilayer, we conducted 60 sets of umbrella-sampling simulations while systematically varying the initial conformation. In each umbrella sampling simulation, each umbrella was simulated for 1.5 μs, yielding a total simulation time of 1.4 ms and 60 independent free energy or potential of mean force (PMF) profiles from bulk water to the center of a POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylcholine) lipid bilayer.The PMF profiles indicate that indolicidin strongly binds to the bilayer, partitioning inside the lipid headgroups (Fig. 1, A and E). Importantly, the mean estimate of ΔGbind0 decays exponentially with equilibration time t_eq, indicating that systematic sampling errors in individual simulations continued to decrease throughout the 1.5-μs interval as rare events led to more favorable states (Fig. 1 B). The low frequency of transitions to more favorable states exacerbates the requirement for massive sampling using multiple independent simulations.Open in a separate windowFigure 1PMF for indolicidin partitioning into a POPC bilayer. (A) Average PMF from 60 independent umbrella-sampling simulations based on 1 < t ≤ 1.5 μs/umbrella. (B) Average ΔGbind0 from 50-ns time intervals per umbrella (t_eq < t ≤ t_eq+50 ns) as a function of equilibration time, t_eq. (Solid line) Single exponential fit to the mean over 0.5 < t_eq ≤ 1.5 μs. (C) Mean values of ΔGbind0 from the 10-μs/umbrella simulations (crosses) together with the mean values of ΔGbind0 (triangles) and exponential fit from panel B. PMF and ΔGbind0 profiles obtained from each of the 60 independent simulations are shown in Fig. S1 in the Supporting Material. (D–F) Representative conformations after 1.5 μs of simulation at zi0 = (D) 3 nm, (E) 1.2 nm, and (F) 0.0 nm. To see this figure in color, go online.Computational limitations precluded extending all 60 sets of umbrella-sampling simulations to even longer times. Instead, we identified the two simulations at each umbrella that appeared to be most representative of equilibrium and extended each to 10 μs per umbrella (see Methods in the Supporting Material). The resulting estimates of ΔGbind0 continued to decrease until t_eq = 4 μs (68 μs in total), after which they stabilized at the asymptotic limit of the exponential fit of the shorter simulation data (ΔGbind0 = −26 ± 5 kcal/mol; Fig. 1 C).As indolicidin approaches the bilayer, it is drawn closer (Fig. 2 A) as salt bridges form between the peptide and the phospholipid headgroups (Fig. 2 B), inducing their protrusion (Figs. 1 D and ​and22 C). At large separation distances, this state is attained only when the peptide becomes highly extended (Fig. 2, D and E). As indolicidin is inserted more deeply, the surface of the lipid bilayer invaginates (Figs. 1 E and ​and22 C), maintaining peptide-lipid salt bridges (Fig. 2 B) and leading to the formation of a pore when the solute is near the bilayer center (Figs. 1 F and ​and22 C, and see Fig. S2, Fig. S3, Fig. S4, and Fig. S5 in the Supporting Material). These Boltzmann-weighted ensemble averages may not be mechanistically representative of nonequilibrium binding events (8,10).Open in a separate windowFigure 2Slow equilibration of bilayer and peptide. (A–D) Color quantifies conformational reorganization for t_eq < t ≤ t_eq + 100 ns as a function of t_eq and |zi0|. (A) Deviation of insertion depth, z, from zi0, Δz ≡ z − zi0; (B) number of peptide-lipid salt bridges, N_SB; (C) volume change of the bilayer’s proximal leaflet in the radial vicinity of the solute, V_ε; and (D) peptide end-to-end distance (EED). There is no sampling for t > 0.5 μs at |zi0| ≥ 4.5 nm. (E) Representative time-series of a trajectory at zi0 = 3.9 nm. (F) Representative conformation at 10 μs for |zi0| = 1.2 nm. To see this figure in color, go online.The reorganization of the peptide, the bilayer, and the ionic interactions between them became more pronounced with increasing simulation time at peptide insertion depths shallower than the global free energy minimum (|zi0| >1.4 nm; Fig. 1 A and Fig. 2, B–D). These conformational transitions are likely the source of the systematic drift of ΔGbind0. Reorganization of the bilayer also controls the rate of equilibration during membrane insertion of an arginine side chain (8,9) and a cyclic arginine nonamer (11), suggesting that slow reorganization of lipids around cationic solutes presents a general impediment to simulation convergence.Consistent with the perturbation of membrane thickness observed by in situ atomic force microscopy (12), our results suggest that indolicidin insertion induces local thinning of the bilayer (Fig. 1, E and F, and Fig. 2, C and F). The different conformational ensembles sampled by the peptide in water and in the lipid bilayer (Fig. 2 D) are consistent with the observations that indolicidin is disordered in solution (13) and adopts stable conformations in the presence of detergent (14). Although the peptide’s conformation continued to change when it was deeply inserted (Fig. 2 D), the amount of water in the bilayer’s hydrophobic core converged relatively rapidly (see Fig. S2). Indolicidin can induce the formation of hydrated, porelike defects (see Fig. S2, Fig. S3, Fig. S4, and Fig. S5) but does not act as a chloride carrier (see Fig. S6 and Fig. S7). Future studies of the mechanism of indolicidin action will examine the effect of multiple peptide binding.The PMF profile presented in this Letter is strikingly different from that computed by Yeh et al. (15) using different force field parameters for indolicidin partitioning into a DMPC (1,2-dimyristoyl-sn-glycero-3-phosphatidylcholine) bilayer, from which the binding free energy was estimated to be 0 kcal/mol (15). However, that study comprised only 25 ns per umbrella and likely suffers from systematic sampling errors induced by initial conditions (see Fig. S8).Our estimate of the binding affinity is much larger than the values obtained for indolicidin and large unilamellar POPC vesicles using isothermal titration calorimetry, −7.4 kcal/mol (16), and equilibrium dialysis, −8.8 kcal/mol (13). Such a discrepancy suggests that the relative accuracy of binding free energies for amino-acid side-chain analogs (8,9,17) does not necessarily extend to polypeptides. Although more work is needed to elucidate the source of this discrepancy, this study underlines the importance of attaining convergence before evaluating force-field accuracy.Importantly, this work also highlights the extensive sampling required to remove systematic errors induced by initial conditions in atomistic simulations of peptides in membranes. Slow equilibration of the system is due to rare transitions across hidden free energy barriers involving reorganization of the membrane. Two simple recommendations are 1), evaluating the time-dependence of ensemble averages, and 2), conducting multiple simulations with different initial conditions. We have recently shown that by using enhanced sampling techniques it is possible to identify the locations of hidden free energy barriers without a priori knowledge (9). Future research will examine strategies for speeding up the crossing of these barriers, such as optimized order parameters including bilayer reorganization and enhanced sampling techniques including a random walk along the order parameter (9).     

DNA,Flexibly Flexible

Investigators have constructed dsDNA molecules with several different base modifications and have characterized their bending and twisting flexibilities using atomic force microscopy, DNA ring closure, and single-molecule force spectroscopy with optical tweezers. The three methods provide persistence length measurements that agree semiquantitatively, and they show that the persistence length is surprisingly similar for all of the modified DNAs. The circular dichroism spectra of modified DNAs differ substantially. Simple explanations based on base stacking strength, polymer charge, or groove occupancy by functional groups cannot explain the results, which will guide further high-resolution theory and experiments.Real double-stranded DNA molecules differ from the idealized zero-Kelvin A, B, and Z forms. They can adopt deformed average conformations, as in bent A-tract DNA or protein-DNA complexes. The path of the DNA helix axis also varies due to thermal energy, so at very long lengths DNA behaves as a random coil. The term “long lengths” is relative to the persistence length P of the wormlike chain model. P is the average offset of the end of a chain along its initial direction, or alternatively the length over which the unit vectors μ1→ and μ2→ tangent to the helix axis lose colinearity according to〈μ1→⋅μ2→〉=〈cosθ〉=e−d12/P,where d₁₂ is the contour length from point 1 to point 2, as in Fig. 1. P can be measured by hydrodynamics (1), atomic force microscopy (AFM) (2), DNA ring closure (3) or protein-DNA looping (4), tethered particle microscopy (5), or single-molecule optical tweezers experiments (6). The long-range loss of memory of DNA direction grows out of local variations in the helix axis direction specified by roll, tilt, and twist angles that parameterize changes in the helix axis direction. For harmonic bending potentials, the bending persistence length is related to roll and tilt according toσroll2+σtilt2=2ℓ/P,where ℓ = 3.4 Å, so for P ∼ 50 nm (147 bp) the average standard deviations in the roll and tilt angles σ_roll and σ_tilt are ∼4.7°, although in real DNA, roll varies more than tilt. Similar relationships hold for twist flexibility (7).Open in a separate windowFigure 1The base modifications studied by Peters et al. (13,14) affect both Watson-Crick hydrogen bonding and groove occupancy. They used AFM, DNA ring closure, single-molecule force spectroscopy, and circular dichroism spectroscopy (not shown) to characterize the resulting changes in bending and twisting flexibility. DNA molecules are not shown to scale. To see this figure in color, go online.DNA flexibility can be studied at contour length scales from Ångstroms to microns. Flexibility at the atomic scale accessed by nuclear magnetic resonance, x-ray crystallography, cryo-electron microscopy, and molecular dynamics simulations (8) refers to many aspects of conformational variability. One active thread of research at this scale concerns interconversion among helical forms, base flipping, DNA kinking, changes in backbone torsion angles, and the sequence dependence of all of these local properties. Local fluctuations in the basepair roll, tilt, and twist angles do seem to predict the correct long-range behavior (9). A second thread asks whether the wormlike chain model holds at DNA lengths shorter than P (2,10); the active controversy concerning enhanced bendability at short lengths has recently been reviewed by Vologodskii and Frank-Kamenetskii (11). A third thread asks whether we can understand the underlying biophysical causes of long-range DNA flexibility. These presumably include base stacking, electrostatic repulsion along the backbone, changes in the counterion atmosphere (12), occupancy of the major and minor grooves by functional groups, conformational entropy, the strength of Watson-Crick hydrogen bonding, and water structure. Helical polymorphisms and the junctions between polymorphs presumably affect the sequence dependence of the persistence length.Peters et al. (13,14) have attempted to understand bending and twisting flexibility by characterizing a variety of modified nucleic acids using DNA ring closure, AFM, and optical tweezer methods, sketched in Fig. 1. In previous work (13), they used ring closure to show that major groove substituents that alter the charge on the polymer do not have substantial effects on the bending persistence length, and that the effects were not correlated in an obvious way to the stacking propensity of the modified bases. The work described in this issue of the Biophysical Journal (14) uses all three methods to demonstrate that DNA with 2-amino-adenosine (a.k.a., 2,6-diaminopurine) substituted for adenosine has an increased persistence length, whereas inosine substitution for guanosine reduces the persistence length, as would be expected if groove occupancy (or the number of Watson-Crick hydrogen bonds) affects flexibility. However, the authors did one experiment too many—when they measured the effects of the earlier major groove substituents (13) using AFM, the correlation with groove occupancy disappeared. This could be because changes in helical geometry, as evidenced by the circular dichroism spectroscopy also reported in the article, alter the grooves sufficiently to prevent a straightforward connection to flexibility.The magnitude of the effect of base modifications on P is the largest for the optical tweezers and the smallest for DNA ring closure, showing that no more than one of the experiments is perfect. The Supporting Material for both articles (13,14) offers valuable resources for the careful evaluation of experimental results and possible sources of error within and between experiments. For example, the DNA lengths and the ionic conditions required by the different methods differ. Ring closure results depend critically on the purity of the DNA and appropriate ligation conditions. Analysis of AFM results averaged several different statistical measures of decaying angular correlations and end-to-end distance, which did not individually always agree. In force spectroscopy there are variations in the bead attachment for each molecule, errors in the stretch modulus can affect the measured persistence length, force can induce DNA melting, and very few molecules can be observed. Rare kinking events proposed to explain enhanced bendability should affect the cyclization experiment most markedly; no evidence for enhanced flexibility was seen. Finally, Peters et al. (14) have observed that DNA twist and twisting flexibility seem to be more sensitive than the persistence length to base modifications.Taken as a whole, this extremely thorough series of experiments shows that we still do not understand the fundamental origins of the remarkable stiffness of double-stranded DNA. There may be compensating effects that make the dissection difficult. For example, changing the charge on the polymer may induce a corresponding adjustment in the counterion condensation atmosphere, leading to a relatively constant residual charge. Groove substituents that enhance basepair stability could enhance bendability for steric reasons. Stacking thermodynamics may not change very much for the very small bend angles at any individual basepair. Locally stiff regions may introduce nearby junctions that are flexible.The stiffness of DNA relative to other biopolymers inspired the development of DNA nanotechnology (although that field has adopted bridged synthetic constructs that are even more rigid). Further research on the biophysics, and specifically the long-range mechanical properties of DNA, will be essential as we build better models of DNA in the cell, which has evolved many proteins that act to increase apparent flexibility. The various aspects of DNA flexibility influence the protein-DNA complexes that mediate DNA’s informational role, the induction of and responses to supercoiling used for long-range communication among sites (15), and chromosome structure and genome organization.     

The gating charge should not be estimated by fitting a two-state model to a Q-V curve

The voltage dependence of charges in voltage-sensitive proteins, typically displayed as charge versus voltage (Q-V) curves, is often quantified by fitting it to a simple two-state Boltzmann function. This procedure overlooks the fact that the fitted parameters, including the total charge, may be incorrect if the charge is moving in multiple steps. We present here the derivation of a general formulation for Q-V curves from multistate sequential models, including the case of infinite number of states. We demonstrate that the commonly used method to estimate the charge per molecule using a simple Boltzmann fit is not only inadequate, but in most cases, it underestimates the moving charge times the fraction of the field.Many ion channels, transporters, enzymes, receptors, and pumps are voltage dependent. This voltage dependence is the result of voltage-induced translocation of intrinsic charges that, in some way, affects the conformation of the molecule. The movement of such charges is manifested as a current that can be recorded under voltage clamp. The best-known examples of these currents are “gating” currents in voltage-gated channels and “sensing” currents in voltage-sensitive phosphatases. The time integral of the gating or sensing current as a function of voltage (V) is the displaced charge Q(V), normally called the Q-V curve.It is important to estimate how much is the total amount of net charge per molecule (Q_max) that relocates within the electric field because it determines whether a small or a large change in voltage is necessary to affect the function of the protein. Most importantly, knowing Q_max is critical if one wishes to correlate charge movement with structural changes in the protein. The charge is the time integral of the current, and it corresponds to the product of the actual moving charge times the fraction of the field it traverses. Therefore, correlating charge movement with structure requires knowledge of where the charged groups are located and the electric field profile. In recent papers by Chowdhury and Chanda (2012) and Sigg (2013), it was demonstrated that the total energy of activating the voltage sensor is equal to Q_max V_M, where V_M is the median voltage of charge transfer, a value that is only equal to the half-point of activation V_1/2 for symmetrical Q-V curves. V_M is easily estimated from the Q-V curve, but Q_max must be obtained with other methods because, as we will show here, it is not directly derived from the Q-V curve in the general case.The typical methods used to estimate charge per molecule Q_max include measurements of limiting slope (Almers, 1978) and the ratio of total charge divided by the number of molecules (Schoppa et al., 1992). The discussion on implementation, accuracy, and reliability of these methodologies has been addressed many times in the literature, and it will not be discussed here (see Sigg and Bezanilla, 1997). However, it is worth mentioning that these approaches tend to be technically demanding, thus driving researchers to seek alternative avenues toward estimating the total charge per molecule. Particularly, we will discuss here the use of a two-state Boltzmann distribution for this purpose. Our intention is to demonstrate that this commonly used method to estimate the charge per molecule is generally incorrect and likely to give a lower bound of the moving charge times the fraction of the field.The two-state Boltzmann distribution describes a charged particle that can only be in one of two positions or states that we could call S₁ and S₂. When the particle with charge Q_max (in units of electronic charge) moves from S₁ to S₂, or vice versa, it does it in a single step. The average charge found in position S₂, Q(V), will depend on the energy difference between S₁ and S₂, and the charge of the particle. The equation that describes Q(V) is:Q(V)=Qmax1+exp[−Qmax(V−V1/2)kT],(1)where V_1/2 is the potential at which the charge is equally distributed between S₁ and S₂, and k and T are the Boltzmann constant and absolute temperature, respectively. The Q(V) is typically normalized by dividing Eq. 1 by the total charge Q_max. The resulting function is frequently called a “single Boltzmann” in the literature and is used to fit normalized, experimentally obtained Q-V curves. The fit yields an apparent V_1/2 (V∼1/2) and an apparent Q_MAX (Q−max), and this last value is then attributed to be the total charge moving Q_max. Indeed, this is correct but only for the case of a charge moving between two positions in a single step. However, the value of Q−max thus obtained does not represent the charge per molecule for the more general (and frequent) case when the charge moves in more than one step.To demonstrate the above statement and also estimate the possible error in using the fitted Q−max from Eq. 1, let us consider the case when the gating charge moves in a series of n steps between n + 1 states, each step with a fractional charge z_i (in units of electronic charge e₀) that will add up to the total charge Q_max.S1↔μ1S2↔μ2…Si↔μiSi+1…Sn↔μnSn+1The probability of being in each of the states S_i is labeled as P_i, and the equilibrium constant of each step is given byμi=exp[zi(V−Vi)kT], i=1…n,where z_i is the charge (in units of e₀) of step i, and V_i is the membrane potential that makes the equilibrium constant equal 1. In steady state, the solution of P_i can be obtained by combiningPi+1Pi=μi, i=1…n and ∑i=1i=n+1Pi=1,givingPi+1=∏m=1iμm1+∑j=1n∏k=1jμk, i=1…n and P1=11+∑j=1n∏k=1jμk.We define the reaction coordinate along the moved charged q asqi=∑j=1izj, i=1…n.The Q-V curve is defined asQ(V)=∑i=1nqiPi+1.Then, replacing P_i yieldsQ(V)=∑i=1n[∑j=1izj][∏m=1iμm]1+∑j=1n∏k=1jμk,or written explicitly as a function of V:Q(V)=∑i=1n[∑j=1izj][∏m=1iexp[zm(V−Vm)kT]]1+∑j=1n∏k=1jexp[zk(V−Vk)kT].(2)Eq. 2 is a general solution of a sequential model with n + 1 states with arbitrary valences and V_i’s for each transition. We can easily see that Eq. 2 has a very different form than Eq. 1, except when there is only a single transition (n = 1). In this latter case, Eq. 2 reduces to Eq. 1 because z₁ and V₁ are equal to Q_max and V_1/2, respectively. For the more general situation where n > 1, if one fits the Q(V) relation obeying Eq. 2 with Eq. 1, the fitted Q−max value will not correspond to the sum of the z_i values (see examples below and Fig. 1). A simple way to visualize the discrepancy between the predicted value of Eqs. 1 and 2 is to compute the maximum slope of the Q-V curve. This can be done analytically assuming that V_i = V_o for all transitions and that the total charge Q_max is evenly divided among those transitions. The limit of the first derivative of the Q(V) with respect to V evaluated at V = V_o is given by this equation:dQ(V)dV|V=V0=Qmax(n+2)12nkT.(3)From Eq. 3, it can be seen that the slope of the Q-V curve decreases with the number of transitions being maximum and equal to Q_max /(4kT) when n = 1 (two states) and a minimum equal to Q_max /(12kT) when n goes to infinity, which is the continuous case (see next paragraph).Open in a separate windowFigure 1.Examples of normalized Q-V curves for a Q_max = 4 computed with Eq. 2 for the cases of one, two, three, four, and six transitions and the continuous case using Eq. 5 (squares). All the Q-V curves were fitted with Eq. 1 (lines). The insets show the fitted valence (Q−max) and half-point (V∼1/2).

Infinite number of steps

Eq. 2 can be generalized to the case where the charge moves continuously, corresponding to an infinite number of steps. If we makez_i = Q_max/n, ?i = 1…n, ??V_i = V_o, ?i = 1…n, then all µ_i = µ, and we can write Eq. 2 as the normalized Q(V) in the limit when n goes to infinity:Qnor(V)=limn→∞∑i=1n[∑j=1iQmaxn]∏m=1iexp[Qmax(V−Vo)nkT]Qmax[1+∑i=1n∏j=1iexp[Qmax(V−Vo)nkT]]=[Qmax(V−Vo)−kT]exp[Qmax(V−Vo)kT]+kTQmax(V−Vo)[exp[Qmax(V−Vo)kT]−1].(4)Eq. 4 can also be written asQnor(V)=12[1+coth[Qmax(V−Vo)2kT]−2kTQmax(V−V0)],(5)which is of the same form of the classical equation of paramagnetism (see Kittel, 2005).

Examples

We will illustrate now that data generated by Eq. 2 can be fitted quite well by Eq. 1, thus leading to an incorrect estimate of the total charge moved. Typically, the experimental value of the charge plotted is normalized to its maximum because there is no knowledge of the absolute amount of charge per molecule and the number of molecules. The normalized Q-V curve, Q_nor, is obtained by dividing Q(V) by the sum of all the partial charges.Fig. 1 shows Q_nor computed using Eq. 2 for one, two, three, four, and six transitions and for the continuous case using Eq. 5 (squares) with superimposed fits to a two-state Boltzmann distribution (Eq. 1, lines). The computations were done with equal charge in each step (for a total charge Q_max = 4e₀) and also the same V_i = −25 mV value for all the steps. It is clear that fits are quite acceptable for cases up to four transitions, but the fit significantly deviates in the continuous case.Considering that experimental data normally have significant scatter, it is then quite likely that the experimenter will accept the single-transition fit even for cases where there are six or more transitions (see Fig. 1). In general, the case up to four transitions will look as a very good fit, and the fitted Q−max value may be inaccurately taken and the total charge transported might be underestimated. To illustrate how bad the estimate can be for these cases, we have included as insets the fitted value of Q−max for the cases presented in Fig. 1. It is clear that the estimated value can be as low as a fourth of the real total charge. The estimated value of V∼1/2 is very close to the correct value for all cases, but we have only considered cases in which all V_i’s are the same.It should be noted that if µ_i of the rightmost transition is heavily biased to the last state (V_i is very negative), then the Q−max estimated by fitting a two-state model is much closer to the total gating charge. In a three-state model, it can be shown that the fitted value is exact when V₁→∞ and V₂→−∞ because in that case, it converts into a two-state model. Although these values of V are unrealistic, the fitted value of Q−max can be very close to the total charge when V₂ is much more negative than V₁ (that is, V₁ >> V₂). On the other hand, If V₁ << V₂, the Q-V curve will exhibit a plateau region and, as the difference between V₁ and V₂ decreases, the plateau becomes less obvious and the curve looks monotonic. These cases have been discussed in detail for the two-transition model in Lacroix et al. (2012).We conclude that it is not possible to estimate unequivocally the gating charge per sensor from a “single-Boltzmann” fit to a Q-V curve of a charge moving in multiple transitions. The estimated Q−max value will be a low estimate of the gating charge Q_max, except in the case of the two-state model or the case of a heavily biased late step, which are rare occurrences. It is then safer to call “apparent gating charge” the fitted Q−max value of the single-Boltzmann fit.

Addendum

The most general case in which transitions between states include loops, branches, and steps can be derived directly from the partition function and follows the general thermodynamic treatment by Sigg and Bezanilla (1997), Chowdhury and Chanda (2012), and Sigg (2013). The reaction coordinate is the charge moving in the general case where it evolves from q = 0 to q = Q_max by means of steps, loops, or branches. In that case, the partition function is given byZ=∑iexp(qi(V−Vi)kT).(6)We can compute the mean gating charge, also called the Q-V curve, asQ(V)=〈q〉=kTZ’Z=kTdlnZdV=∑iqiexp(qi(V−Vi)kT)∑iexp(qi(V−Vi)kT).(7)The slope of the Q-V is obtained by taking the derivative of 〈q〉 with respect to V:dQ(V)dV=(kT)2d2lnZdV2.(8)Let us now consider the gating charge fluctuation. The charge fluctuation will depend on the number of possible conformations of the charge and is expected to be a maximum when there are only two possible charged states to dwell. As the number of intermediate states increases, the charge fluctuation decreases. Now, a measure of the charge fluctuation is given by the variance of the gating charge, which can be computed from the partition function as:〈Δq2〉=〈q2〉−〈q〉2=(kT)2(Z’’Z−(Z’Z)2)=(kT)2d2lnZdV2.(9)But the variance (Eq. 9) is identical to the slope of Q(V) (Eq. 8). This implies that the slope of the Q-V is maximum when there are only two states.     

The limitations,dangers, and benefits of simple methods for testing identifiability

In their Commentary paper, Villaverde and Massonis (On testing structural identifiability by a simple scaling method: relying on scaling symmetries can be misleading) have commented on our paper in which we proposed a simple scaling method to test structural identifiability. Our scaling invariance method (SIM) tests for scaling symmetries only, and Villaverde and Massonis correctly show the SIM may fail to detect identifiability problems when a model has other types of symmetries. We agree with the limitations raised by these authors but, also, we emphasize that the method is still valuable for its applicability to a wide variety of models, its simplicity, and even as a tool to introduce the problem of identifiability to investigators with little training in mathematics.

In their Commentary paper, Villaverde and Massonis (On testing structural identifiability by a simple scaling method: relying on scaling symmetries can be misleading [1]) have commented on our paper in which we proposed a simple scaling method to test structural identifiability [2]. Our scaling invariance method (SIM) tests for scaling symmetries only, and Villaverde and Massonis correctly show the SIM may fail to detect identifiability problems when a model has other types of symmetries (we indeed indicated but not investigated the importance of generalizing the method to other symmetries). Thus, we agree that our simple method provides a necessary but not sufficient condition for identifiability, and we appreciate their careful analysis and constructive criticism.We nevertheless think that the simple method remains useful because it is so simple. Even for investigators with little training in mathematics, the method provides a necessary condition for structural identifiability that can be derived in a few minutes with pen and paper. Similarly, we have found its pedagogic strength by teaching the method to our own graduate students and colleagues. More advanced methods (such as STRIKE-GOLDD [3,4], COMBOS [5], or SIAN [6]) are typically intimidating for researchers with a background in Biology or Bioinformatics. This simple method can help those practitioners to familiarize themselves with the identifiability problem and better understand their models.Finally, it is worth noting that if scaling invariance is the only symmetry (as it was in all the cases we analyzed), our SIM remains valuable (albeit uncontrolled), and surprisingly effective for a wide variety of problems (as the extensive list collected in the Supplementary Material our paper [2]). We guess that the SIM especially fails when applied to linear models (as more potential rotations of the variables leave the system invariant), and in non-linear scenarios where some parameters are identical. For instance, the FitzHugh-Nagumo model raised by Villaverde and Massonis, x˙1(t)=c(x1(t)−x13(t)3−x2(t)+d),x˙2(t)=1c(x1(t)+a−b·x2(t)),y(t)=x1(t), could have been written as x˙1(t)=λ1x1(t)−λ2x13(t)3−λ3x2(t)+d,x˙2(t)=λ4x1(t)+a−b·x2(t),y(t)=x1(t) where λ₁ = λ₂ = λ₃ = 1/λ₄ = c. One of the reasons why our method fails, in this case, might be these additional symmetries introduced in this more elaborate notation of the model.Hence, it is worth understanding generic conditions under which the SIM method is expected to be fragile, possibly using STRIKE-GOLDD to test large families of nonlinear models.As a final remark, we appreciate that Villaverde and Massonis have shared their source code, so researchers might have a gold standard to test identifiability.     

Detecting Subtle Plasma Membrane Perturbation in Living Cells Using Second Harmonic Generation Imaging

The requirement of center asymmetry for the creation of second harmonic generation (SHG) signals makes it an attractive technique for visualizing changes in interfacial layers such as the plasma membrane of biological cells. In this article, we explore the use of lipophilic SHG probes to detect minute perturbations in the plasma membrane. Three candidate probes, Di-4-ANEPPDHQ (Di-4), FM4-64, and all-trans-retinol, were evaluated for SHG effectiveness in Jurkat cells. Di-4 proved superior with both strong SHG signal and limited bleaching artifacts. To test whether rapid changes in membrane symmetry could be detected using SHG, we exposed cells to nanosecond-pulsed electric fields, which are believed to cause formation of nanopores in the plasma membrane. Upon nanosecond-pulsed electric fields exposure, we observed an instantaneous drop of ∼50% in SHG signal from the anodic pole of the cell. When compared to the simultaneously acquired fluorescence signals, it appears that the signal change was not due to the probe diffusing out of the membrane or changes in membrane potential or fluidity. We hypothesize that this loss in SHG signal is due to disruption in the interfacial nature of the membrane. The results show that SHG imaging has great potential as a tool for measuring rapid and subtle plasma membrane disturbance in living cells.As the epicenter for many cellular functions, understanding the dynamics of the plasma membrane is important to monitoring biological phenomena. External forces acting upon the plasma membrane (e.g., electric, mechanical) have been shown to cause rapid disturbances, often resulting in dramatic changes in cell physiology (1–3). To understand this interaction, a minimally invasive, highly sensitive imaging technique that enables monitoring the structure of the plasma membrane is needed. Lipophilic dyes, which embed themselves into lipid membranes, are sensitive to the surrounding electric field and, therefore, report changes in membrane fluidity as well as voltage due to the capacitive nature of the membranes (4,5). This sensitivity is typically detected as a shift in the fluorescence emission spectrum. Localization of the fluorescence signal to only the plasma membrane is difficult because the probes also label internal membrane structures. Thus, to overcome this lack of spatial selectivity, second harmonic generation (SHG) has been used as an alternative to fluorescence for membrane imaging (6,7).In SHG, a second-order nonlinear polarization is induced by electronic disruption of a probe molecule from the electromagnetic field of the incident laser beam. This polarization generates oscillating dipole moments that reradiate light at twice the energy of the excitation beam. The induction of this dipole is sensitive to the static electric field surrounding the probe and the steady-state molecular polarization of the probe molecule. These properties make SHG probes useful for monitoring changes in biological membranes.First, as the voltage potential across the membrane changes, the static electric field around the probe also shifts, making the probe sensitive to these variations (7). Several SHG probes have, therefore, been employed to monitor plasma membrane potential (7,8).Second, because the dipole is affected by the steady-state molecular polarization of the probe itself, a SHG signal is only produced in materials that lack a center of inversion symmetry. In the centrosymmetric case, any emitted radiation is cancelled out by destructive interference. The properties of an interfacial environment, such as a cellular plasma membrane, not only provide the necessary asymmetry, but cause the polarized lipophilic dyes to be aligned in respect to the interface, instead of being randomly distributed as they would in a bulk environment. This alignment allows the generation of a coherent SHG signal from the plasma membrane while the rest of the cell remains nearly signal-free (6,7).We investigated whether the alignment sensitivity of the SHG response could be used to detect minute changes in the organization of the plasma membrane. Jurkat clone E6-1 human T-lymphocytes with a spherical morphology were selected for optimum signal clarity and cultured as directed by American Type Culture Collection (ATCC, Manassas, VA) with 1 I.U./mL penicillin and 0.1 μg/mL streptomycin. Cells were added to 35-mm poly-L-lysine-coated glass-bottomed dishes (MatTek, Ashland, MA) and incubated for 1 h in growth media to allow adherence. Before loading, the cells were rinsed with a buffer consisting of 135 mM NaCl, 5 mM KCl, 2 mM MgCl2, 10 mM HEPES, 10 mM glucose, 2 mM CaCl2, pH 7.4, 290–310 mOsm. SHG probes, Di-4-ANEPPDHQ (Di-4) (5 μM final concentration), FM4-64 (15 μM) or ATR (100 μM, 1 mg/mL BSA) were added to the buffer solution and incubated for 1 h. Cellular imaging was performed in the labeling buffer to limit diffusion of the probe molecules out of the cell membranes.A Ti:sapphire oscillator at 980 nm (Coherent Chameleon, 130 fs, 80 MHz, ∼15 mW at the sample; Coherent Laser, Santa Clara, CA) was coupled through the scan head of a modified model No. TCS SP5 II (Leica Geosystems, Norcross, GA) for SHG and multiphoton-excited fluorescence imaging (40×, water, 1.1 NA) using resonant scanning. SHG signal was collected in transmission by a photomultiplier tube after 680-nm shortpass and 485/25-nm bandpass filters; simultaneous fluorescence signal was collected in the epi-direction by two non-descanned photomultiplier tubes with 540/60-nm and 650/60-nm bandpass filters.Three SHG probes previously used to monitor voltage or membrane order in living cells were tested (8–10). Although ATR is reported to be effective in monitoring membrane voltage, we obtained nearly no SHG signal, despite successful loading as indicated by the fluorescence signal (Fig. 1). When FM4-64 and Di-4 were loaded to similar fluorescence intensities, nearly equivalent SHG signal was collected. Di-4 did appear to have a greater internalization of the dye. However, after the first frame, the FM4-64 signal dropped considerably (Fig. 1 b), an observation reported as a membrane voltage-independent bleaching effect (8). This drop in signal recovered after excitation was blocked for several seconds, but quantification of the response was difficult. Di-4 did not suffer as dramatic a drop in signal upon excitation, and still had sufficient SHG signal/noise after several seconds, so it was used in all further experiments.Open in a separate windowFigure 1(a) Fluorescence (top) and SHG (bottom) images for the three probes. (b) Signal/noise for the fluorescence and SHG for the initial frame and shortly after beginning acquisition. Error bars represent the mean ± SE (n = 10). Scale bar is 10 μm.To test whether Di-4 would report a rapid change in membrane organization, we applied a single nanosecond-duration pulsed electric field (nsPEF) to the labeled cell. These ultrabrief, high-intensity (MV/m) pulses differ from longer (μs-ms), lower-intensity (kV/m) pulses traditionally associated with electroporation in induced cellular response (3,11,12). Through selective uptake of small ions (Ca²⁺, Ti⁺) with limited uptake of propidium iodide, nsPEF have been previously postulated to cause nanopores (<2 nm diameter) in the plasma membrane. In contrast with a previous study observing poration resulting from traditional electroporation (13), the brevity of this apparently novel cellular insult allows for the decoupling of the mechanical effects of the pulse on the membrane from the electrical effects of the pulse itself. A single pulse, generated by a custom pulse generator, was delivered to the cells using a pair of 125-μm diameter tungsten electrodes, separated edge-to-edge by 150 μm, as previously described in Ibey et al. (14). For maximum visualization of changes in the SHG signal, a half-wave plate was placed before the scan head to align the polarization of the laser such that the brightest signals from the plasma membrane were at the poles facing the electrodes.The Di-4 SHG signal in response to a single 16.7 kV/cm, 600-ns nsPEF is shown in Fig. 2. Before the pulse, the intensity of the SHG signal is high at each of these poles. Immediately after the pulse, the SHG intensity drops by ∼50% on the side of the cell facing the anodic electrode, whereas little intensity is lost at the other pole. This response is plotted in Fig. 2 (pulse applied at 2 s), where it is apparent that the response is near instantaneous with little recovery in signal in the 5 s postexposure. The SHG response matches the previously observed effect of this stimulus, where ion uptake displayed a polar dependence and persisted for a number of minutes (11,12). Images taken 5 min after an nsPEF exposure are also shown in Fig. 2. These images confirm the eventual recovery of the cell and the corresponding return of SHG signal to preexposure levels.Open in a separate windowFigure 2(a) SHG images showing drop in signal on the anodic (or A-pole) of the cell. (b) Time trace of SHG response with the electrical pulse applied at 2 s that shows a near-instantaneous drop in the SHG signal at the anodic pole of the cell. (c) SHG image preexposure, immediately postexposure, and then 5-min postexposure showing recovery of the SHG signal.To decouple membrane disturbance from environmental changes around the membrane, we compared the SHG response to the simultaneously acquired fluorescence signal. Because fluorescence is not subject to the strict orientation requirement of SHG, the plasma membrane fluorescence signal provides an indication of the membrane fluidity and/or potential. Despite the dramatic shift in SHG intensity on the anodic pole upon the electrical pulse exposure (Fig. 3 a), the fluorescence channels display little response from the equivalent membrane sections with the exception of photobleaching and a slight increase in signal in both emission bands on the anodic side (Fig. 3, b and c). The shading in these graphs represents the mean ± SE for six cells. Although this slight increase may indicate that a small amount of dye is simply diffusing in or out of the membrane upon exposure, the fluorescence response is not as rapid or as lasting as the SHG response. Change in membrane fluidity or voltage can also be quantized using these fluorescence signals and a value known as the generalized polarization (GP) (4),GP=I515−570−I620−680I515−570+I620−680.(1)As with the raw intensity of the individual signals, the GP value for the membrane (Fig. 3 d) shows no significant shift, indicating that the membrane is likely not transitioning between a more raft- and fluidlike state. Thus, it seems likely that the dye was initially aligned in the tightly-packed ordered membrane so that the probes were able to generate a SHG photon. As shown in Fig. 3 e, we postulate that upon electrical pulse exposure, the membrane was disrupted by the formation with nanopores giving the probe molecules the flexibility to disorient within the membrane. The resulting alignment of the probes is more isotropic in nature, thereby limiting the probes probability of producing a SHG photon. The fluorescence signal remained, however, indicating that the probes remained active in the membrane.Open in a separate windowFigure 3(a) Average SHG signal showing the dramatic drop in signal on the anodic pole at the pulse application (2 s). (b and c) Simultaneous TPF signals showing nearly no instantaneous change at the pulse application. (d) GP showing no observable changes in the membrane potential or fluidity after the pulse. (Shaded areas) Fit to the mean ± SE for each trace (n = 6). (e) Conceptualization of the hypothesized membrane disruption underlying the observed change in SHG response.Thus, by taking advantage of the selection criteria of SHG, we were able to successfully use the SHG probe, Di-4, to monitor rapid disruption of the plasma membrane. Because SHG can only be generated when the probes are aligned in the plasma membrane, the SHG signal diminishes significantly upon disruption. The simultaneous collection of the multiphoton-excited fluorescence signal was advantageous in that it demonstrated that the probes did not simply diffuse out of the membrane, did not appear to be energetically disrupted by the electric pulse, and showed that the membrane changes were not simply a change in lipid order. We believe that this technique holds tremendous potential for use in the study of how external stimuli interact with and change the orientation of biological membranes. Such knowledge may allow for further understanding of how manipulation of cells and biological systems can be achieved using external stimuli.     

Type-3 BRET,an Improved Competition-Based Bioluminescence Resonance Energy Transfer Assay

The heart adjusts its power output to meet specific physiological needs through the coordination of several mechanisms, including force-induced changes in contractility of the molecular motor, the β-cardiac myosin (βCM). Despite its importance in driving and regulating cardiac power output, the effect of force on the contractility of a single βCM has not been measured. Using single molecule optical-trapping techniques, we found that βCM has a two-step working stroke. Forces that resist the power stroke slow the myosin-driven contraction by slowing the rate of ADP release, which is the kinetic step that limits fiber shortening. The kinetic properties of βCM are affected by load, suggesting that the properties of myosin contribute to the force-velocity relationship in intact muscle and play an important role in the regulation of cardiac power output.The cardiac cycle is a tightly regulated process in which the heart generates power during systole and relaxes during diastole. Appropriate power must be generated to effectively pump blood against cardiac afterload. Dysfunction of this cycle has devastating consequences for affected individuals.Cardiac power output is regulated by several feedback mechanisms (e.g., neuronal, hormonal, mechanical) that ultimately lead to changes in the force and power output of the molecular motor, β-cardiac myosin (βCM). In isolated cardiac fibers and cardiomyocytes, loading the muscle during systole slows contraction and alters power output. It is widely believed that this slowing is partially due to force-induced inhibition of myosin ATPase kinetics, similar to the Fenn Effect in skeletal muscle. However, this hypothesis has not been directly tested at the molecular level. Much of our contemporary view of how power is generated in cardiac muscle is due to in vivo and isolated muscle-fiber studies (1). Substantial progress has been made in understanding the actomyosin interactions required for power generation, but resolving the molecular effects of mechanical load on the ATPase properties of βCM in intact muscle has been challenging. Nevertheless, determining the biophysical parameters that define βCM contractility is key to understanding cardiac regulation and the etiology of several muscle diseases (1).In vitro assays using isolated contractile proteins have been central to advancing our understanding of contractility, although most experiments have been conducted at low resisting loads that do not mimic working conditions. Elegant optical trapping experiments have imposed loads on small ensembles of murine α-cardiac myosin at subsaturating [ATP] (2), and these experiments suggest that force slows α-cardiac myosin kinetics. The kinetic properties of α-cardiac myosin are substantially different from βCM, the primary isoform in the adult human myocardium (3). Thus, experiments using βCM must be performed to determine the unitary force-dependent kinetic parameters of this key molecular motor. We used optical trapping to measure the working-stroke displacement and force dependence of actin-detachment kinetics of single porcine βCM molecules at saturating ATP concentrations. These experiments allow direct measurement of the force-velocity relationship for single βCM molecules and reveal the mechanism of how loads regulate βCM-driven power output.Using the three-bead geometry (4) in which an actin filament is strung between two optically-trapped beads and then lowered over a bead that is sparsely coated with purified full-length porcine ventricular βCM, interactions between single βCM molecules and actin were recorded at 10 μM ATP (Fig. 1 A) (5, 6). Ensemble averages of these interactions were constructed to determine the size and kinetics of the working stroke (7, 8). βCM has an average displacement (6.8 ± 0.04 nm) that is similar to previously characterized muscle myosins (9, 10). Similar to skeletal muscle myosin (10), ensemble averages clearly show that the βCM working stroke is composed of two substeps with average displacements of 4.7 ± 0.05 nm and 1.9 ± 0.05 nm (Fig. 1 B). A single exponential function was fit to the rising-phase of the time-forward ensemble averages, yielding a rate (74 ± 2 s⁻¹) for the transition from state 1 to state 2 (Fig. 1 C). This rate is similar to the biochemical rate of ADP release measured for βCM (64 s⁻¹) (3), indicating that this structural transition is associated with the release of ADP. The rate of the rising phase of the time-reversed ensemble averages (22 ± 0.7 s⁻¹) reports the rate of exit from state 2 and is consistent with the biochemical rate of ATP binding and actomyosin detachment at 10 μM ATP (16 s⁻¹) (3) (Fig. 1 C).Open in a separate windowFigure 1(A) Representative data trace showing actomyosin displacements generated by βCM at 10 μM ATP. (Blue lines) Individual binding events. (B) Ensemble averages of the βCM working stroke generated from averaging 1295 binding interactions collected at 10 μM ATP. Single exponential functions were fit to the data (red lines) and the reported errors are the standard errors from the fit. (C) Cartoon showing an idealized actomyosin interaction with the corresponding mechanical and biochemical states. To see this figure in color, go online.To examine actomyosin detachment kinetics under working conditions, a positional feedback optical clamp was used to apply a dynamic load to the myosin, keeping the myosin at an isometric position during its working stroke (11). We measured the effect of force on the actin-attachment duration at 4 mM Mg.ATP to ensure that the rate of ATP binding is not rate-limiting for detachment. Increases in attachment durations are observed as the force on the myosin is increased (Fig. 2 A, inset). Assuming a two-state model (12), we expect the attachment durations to be exponentially distributed at each force with the force-dependent actin detachment rate, k(F), given by (13)k(F)=k0∗e−F·ddetkB∗T,(1)where k₀ is the rate of the primary force-sensitive transition in the absence of force, F is the force on the myosin, d_det is the distance to the transition state (a measurement of force sensitivity), k_B is Boltzmann’s constant, and T is the temperature. Maximum likelihood estimation (MLE) fitting Eq. 1 to the data yields a detachment rate (k₀ = 71 (−1.0/+0.8 s⁻¹)) that is similar to the rate of ADP release measured for βCM (64 s⁻¹) (3) and the rate of the time-forward ensemble averages (74 ± 2 s⁻¹). Thus, the ADP release step (and the accompanying state-1 to state-2 mechanical transition) is force-sensitive (d_det = 0.97 (−0.014/+0.011) nm). The value of d_det indicates that the ADP release step slows with increasing force, but less than some other characterized myosins (14). Using the values determined from the MLE fitting and the measured size of the working stroke, it is possible to calculate a force-velocity relationship for βCM, assuming the rate of ADP release limits actin motility (Fig. 2 A).Open in a separate windowFigure 2(A, Inset) Single molecule actomyosin interactions were collected in the presence of the isometric optical clamp. The scatter plot shows 262 binding events. Attachment durations are exponentially distributed at each force. (A) The detachment rate as a function of force as determined by MLE fitting. (Black line) Best fit; (small gray shaded area) 95% confidence interval. (Right axis) Velocity, calculated by multiplying the displacement of the working stroke by the detachment rate. (B) The calculated mean detachment rate as a function of force. Attachment durations were binned according to the average force experienced by the myosin during the binding event. Error bars were calculated via bootstrapping simulations of each force bin. (Blue line) Expected mean detachment rate based on the MLE fitting and the limited temporal resolution of our experiment (see the Supporting Material for details). (C) Proposed model for how force slows shortening velocity. Force inhibits the mechanical transition associated with ADP release, slowing the rate of actomyosin detachment. To see this figure in color, go online.The MLE fitting of Eq. 1 assumes an exponential distribution of attachment durations at every force. As such, the MLE fitting of the raw data should yield correct values of the parameters k₀ and d_det, despite limitations of the temporal resolution of our experiment (see Supporting Material for detailed discussion of MLE fitting). Frequently, groups report the mean attachment duration as a function of force. However, the mean attachment duration at each force will be overestimated because some shorter binding events cannot be resolved. We provide a method for calculating the expected mean detachment rate based on the parameters determined from the MLE fitting, given the limited temporal resolution of the experiment, and verify the robustness of the MLE fitting (see the Supporting Material). For demonstration purposes only, Fig. 2 B shows that the measured mean detachment rate agrees well with the expected mean detachment rate based on the MLE fitting and the temporal resolution of the experiment. It should be emphasized that the relevant dissociation values are obtained from the MLE fitting in Fig. 2 A (see also Figs. S1–S3).Our data demonstrate that at saturating [ATP], the detachment rate is limited by the ADP release step, which is the same transition that limits fiber shortening velocity (15). We propose that resisting loads slow ADP release and actin detachment by slowing the mechanical transition that accompanies ADP release (Fig. 2 C), thereby reducing the shortening velocity of muscle fibers. Thus, our data demonstrate that the intrinsic force-dependent properties of βCM contribute to the force-velocity relationship in the heart. It is important to note that our proposed mechanism does not rule out additional mechanisms by which force could directly modulate the activity of actomyosin such as force-induced reversal of the power stroke (11) or population of branched pathways (16, 17).Are the loads in our experiments physiologically relevant to contracting muscle? Modeling of the force per cross-bridge generated in isometric soleus muscle, which contains the βCM isoform, suggests a load of 2–4 pN per myosin (18). At these loads, we expect actin-detachment to slow up to threefold. Interestingly, βCM is substantially less force-sensitive than smooth muscle myosin (d_det = 2.7), suggesting that βCM can generate more power (the product of force and velocity) under load.In conclusion, our data show that cardiac power output can be directly modulated by force at the level of single myosin molecules. These data will enable the comparison of how molecular changes, such as light-chain phosphorylation, pharmacological treatments, or mutations associated with cardiomyopathies, affect the ability of the myosin to generate power against the afterload.     

Mutations in PBP2 from ceftriaxone-resistant Neisseria gonorrhoeae alter the dynamics of the β3–β4 loop to favor a low-affinity drug-binding state

Resistance to the extended-spectrum cephalosporin ceftriaxone in the pathogenic bacteria Neisseria gonorrhoeae is conferred by mutations in penicillin-binding protein 2 (PBP2), the lethal target of the antibiotic, but how these mutations exert their effect at the molecular level is unclear. Using solution NMR, X-ray crystallography, and isothermal titration calorimetry, we report that WT PBP2 exchanges dynamically between a low-affinity state with an extended β3–β4 loop conformation and a high-affinity state with an inward β3–β4 loop conformation. Histidine-514, which is located at the boundary of the β4 strand, plays an important role during the exchange between these two conformational states. We also find that mutations present in PBP2 from H041, a ceftriaxone-resistant strain of N. gonorrhoeae, increase resistance to ceftriaxone by destabilizing the inward β3–β4 loop conformation or stabilizing the extended β3–β4 loop conformation to favor the low-affinity drug-binding state. These observations reveal a unique mechanism for ceftriaxone resistance, whereby mutations in PBP2 lower the proportion of target molecules in the high-affinity drug-binding state and thus reduce inhibition at lower drug concentrations.Keywords: PBP2, Neisseria gonorrhoeae, beta-lactam, conformational dynamics, antibiotic resistance

Neisseria gonorrhoeae is the causative agent of the sexually transmitted infection gonorrhea, with nearly 80 million cases worldwide each year (1). Without antibiotic treatment, infections persist as a chronic disease and can cause serious sequelae, including pelvic inflammatory disease, infertility, arthritis, and disseminated infections (2). For many years, N. gonorrhoeae was treated with a single dose of penicillin, and more recently, ceftriaxone. In 2012, the emergence of several high-level ceftriaxone-resistant strains led the Centers for Disease Control and Prevention to change its recommended treatment for gonorrhea from monotherapy to dual therapy with ceftriaxone and azithromycin (3, 4, 5). However, treatment failures have been reported for both agents, and in 2018, a strain with high-level resistance to both ceftriaxone and azithromycin was identified (6, 7). Concern about azithromycin resistance led the Centers for Disease Control and Prevention recently to drop the recommendation of dual therapy in favor of an increased dose (500 mg) of ceftriaxone alone (8). Both penicillin and ceftriaxone inhibit cell wall biosynthesis in N. gonorrhoeae by targeting penicillin-binding protein 2 (PBP2).PBP2 is an essential peptidoglycan transpeptidase (TPase) that crosslinks the peptide chains from adjacent peptidoglycan strands during cell-wall synthesis (9). β-lactam antibiotics, including the extended-spectrum cephalosporin (ESC) ceftriaxone, are analogs of the d-Ala-d-Ala C terminus of the peptidoglycan substrate and as such target PBP2 by binding to and reacting with the active-site serine nucleophile (Ser310 in N. gonorrhoeae PBP2) to form a covalently acylated complex (10, 11). The acylation reaction (Equation 1) proceeds first through formation of a noncovalent complex with the β-lactam (defined by the equilibrium constant, Ks), which is then attacked by the serine nucleophile to form a covalent acyl-enzyme complex (k₂). For PBPs, hydrolysis of the acyl-enzyme (k₃) is very slow compared with its formation, and the enzyme is essentially irreversibly inactivated. The acylation of PBPs by β-lactam antibiotics is therefore defined by a second-order rate constant, k₂/Ks (M⁻¹ s⁻¹), which reflects both the noncovalent binding affinity (Ks) and the first-order acylation rate (k₂):E+S↔KsE⋅S→k2E−S→k3E+P(1)The emergence of resistance to penicillin and ceftriaxone in N. gonorrhoeae occurs primarily via the acquisition of mutant alleles of the penA gene encoding PBP2 (12). These alleles are referred to as mosaic because they arise through multiple homologous recombination events with DNA released by commensal Neisseria species. PBP2 from the high-level ceftriaxone-resistant strain, H041, contains 61 mutations compared with PBP2 from the antibiotic-susceptible strain, FA19 (13, 14). Determining how these mutations lower the k₂/Ks of ceftriaxone for PBP2 by over 10,000-fold while still preserving essential TPase activity is fundamental for understanding the evolution of antibiotic resistance.Toward this goal, we have identified a subset of these mutations that, when incorporated into the penA gene from FA19, confer ∼80% of the increase in minimum inhibitory concentration for ceftriaxone relative to that of the penA gene from H041 (penA41) (15, 16). We recently reported the structures of apo and ceftriaxone-acylated PBP2 at high resolution and have detailed conformational changes in β3 and the β3–β4 loop involved in antibiotic binding and acylation (17). Intriguingly, although present in the active site region, most of the mutations conferring resistance are not in direct contact with ceftriaxone in the crystal structure of acylated PBP2 (17, 18). We have proposed that these mutations alter the binding and acylation kinetics of PBP2 with ceftriaxone by restricting protein dynamics (18).To understand further the structural and biochemical mechanisms by which these mutations lower the acylation rates of β-lactam antibiotics, we utilized a combination of solution ¹⁹F NMR, X-ray crystallography, and biochemical approaches to investigate PBP2. We report that the β3–β4 loop in the TPase domain of WT PBP2, which is known to adopt markedly different conformations in the apo versus acylated crystal structures (17), samples two major conformational states in solution. Substitutions of WT PBP2 residues with mutations in H041 that confer ceftriaxone resistance alter the conformational landscape of PBP2 by destabilizing the high-affinity state containing the inward conformation of the β3–β4 loop and stabilizing a low-affinity conformation containing an extended β3–β4 loop conformation, thereby restricting access to the inward conformation required for high-affinity drug binding. Our combined solution NMR and crystallographic analyses of PBP2 and its preacylation drug complexes further support the notion that mutations in PBP2 from ceftriaxone-resistant strains of N. gonorrhoeae confer antibiotic resistance by hindering conformational changes required to form a productive drug-binding state (18).