Dissociation of the Octameric Enolase from S. Pyogenes - One Interface Stabilizes Another

Most enolases are homodimers. There are a few that are octamers, with the eight subunits arranged as a tetramer of dimers. These dimers have the same basic fold and same subunit interactions as are found in the dimeric enolases. The dissociation of the octameric enolase from S. pyogenes was examined, using NaClO₄, a weak chaotrope, to perturb the quaternary structure. Dissociation was monitored by sedimentation velocity. NaClO₄ dissociated the octamer into inactive monomers. There was no indication that dissociation of the octamer into monomers proceeded via formation of significant amounts of dimer or any other intermediate species. Two mutations at the dimer-dimer interface, F137L and E363G, were introduced in order to destabilize the octameric structure. The double mutant was more easily dissociated than was the wild type. Dissociation could also be produced by other salts, including tetramethylammonium chloride (TMACl) or by increasing pH. In all cases, no significant amounts of dimers or other intermediates were formed. Weakening one interface in this protein weakened the other interface as well. Although enolases from most organisms are dimers, the dimeric form of the S. pyogenes enzyme appears to be unstable.     

A Model of the Quaternary Structure of Enolases, Based on Structural and Evolutionary Analysis of the Octameric Enolase from Bacillus subtilis

Purified enolase from Bacillus subtilis has a native mass of approximately 370 kDa. Since B. subtilis enolase was found to have a subunit mass of 46.58 kDa, the quaternary structure of B. subtilis is octameric. The pl for B. subtilis enolase is 6.1, the pH optimum (pH_o) for activity is 8.1–8.2, and the K _m for 2-PGA is approximately 0.67 mM. Using the dimeric C structure of yeast dimeric enolase as a guide, these dimers were arranged as a tetramer of dimers to simulate the electron microscopy image processing obtained for the octameric enolase purified from Thermotoga maritima. This arrangement allowed identification of helix J of one dimer (residues 86–96) and the loop between helix L and strand 1 (HL–S1 loop) of another dimer as possible subunit interaction regions. Alignment of available enolase amino acid sequences revealed that in 16 there are two tandem glycines at the C-terminal end of helix L and the HL–S1 loop is truncated by 4–6 residues relative to the yeast polypeptide, two structural features absent in enolases known to be dimers. From these arrangements and alignments it is proposed that the GG tandem at the C-terminal end of helix L and truncation of the HL–S1 loop may play a critical role in octamer formation of enolases. Interestingly, the sequence features associated with dimeric quaternary structure are found in three phylogenetically disparate groups, suggesting that the ancestral enolase was an octamer and that the dimeric structure has arisen independently multiple times through evolutionary history.     

Dominance of Deleterious Alleles Controls the Response to a Population Bottleneck

Population bottlenecks followed by re-expansions have been common throughout history of many populations. The response of alleles under selection to such demographic perturbations has been a subject of great interest in population genetics. On the basis of theoretical analysis and computer simulations, we suggest that this response qualitatively depends on dominance. The number of dominant or additive deleterious alleles per haploid genome is expected to be slightly increased following the bottleneck and re-expansion. In contrast, the number of completely or partially recessive alleles should be sharply reduced. Changes of population size expose differences between recessive and additive selection, potentially providing insight into the prevalence of dominance in natural populations. Specifically, we use a simple statistic, BR≡∑xipop1/∑xjpop2, where x _i represents the derived allele frequency, to compare the number of mutations in different populations, and detail its functional dependence on the strength of selection and the intensity of the population bottleneck. We also provide empirical evidence showing that gene sets associated with autosomal recessive disease in humans may have a B _R indicative of recessive selection. Together, these theoretical predictions and empirical observations show that complex demographic history may facilitate rather than impede inference of parameters of natural selection.     

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.     

On the Treewidths of Graphs of Bounded Degree

In this paper, we develop a new technique to study the treewidth of graphs with bounded degree. We show that the treewidth of a graph G = (V, E) with maximum vertex degree d is at most (1−Ce−4.06d)|V| for sufficiently large d, where C is a constant.     

Enhancement of Chemotactic Cell Aggregation by Haptotactic Cell-To-Cell Interaction

The crawling of biological cell is a complex phenomenon involving various biochemical and mechanical processes. Some of these processes are intrinsic to individual cells, while others pertain to cell-to-cell interactions and to their responses to extrinsically imposed cues. Here, we report an interesting aggregation dynamics of mathematical model cells, when they perform chemotaxis in response to an externally imposed global chemical gradient while they influence each other through a haptotaxis-mediated social interaction, which confers intriguing trail patterns. In the absence of the cell-to-cell interaction, the equilibrium population density profile fits well to that of a simple Keller-Segal population dynamic model, in which a chemotactic current density J→chemo∼∇p competes with a normal diffusive current density J→diff∼∇ρ, where p and ρ refer to the concentration of chemoattractant and population density, respectively. We find that the cell-to-cell interaction confers a far more compact aggregation resulting in a much higher peak equilibrium cell density. The mathematical model system is applicable to many biological systems such as swarming microglia and neutrophils or accumulating ants towards a localized food source.     

Hierarchical Model of Fibrillar Collagen Organization for Interpreting the Second-Order Susceptibility Tensors in Biological Tissue

The second-order nonlinear polarization properties of fibrillar collagen in various rat tissues (vertebrae, tibia, tail tendon, dermis, and cornea) are investigated with polarization-dependent second-harmonic generation (P-SHG) microscopy. Three parameters are extracted: the second-order susceptibility ratio, R = χZZZ(2)′/χZXX(2)′; a measure of the fibril distribution asymmetry, |A|; and the weighted-average fibril orientation, 〈δ〉. A hierarchical organizational model of fibrillar collagen is developed to interpret the second-harmonic generation polarization properties. Highlights of the model include: collagen type (e.g., type-I, type-II), fibril internal structure (e.g., straight, constant-tilt), and fibril architecture (e.g., parallel fibers, intertwined, lamellae). Quantifiable differences in internal structure and architecture of the fibrils are observed. Occurrence histograms of R and |A| distinguished parallel from nonparallel fibril distributions. Parallel distributions possessed low parameter values and variability, whereas nonparallel distributions displayed an increase in values and variability. From the P-SHG parameters of vertebrae tissue, a three-dimensional reconstruction of lamellae of intervertebral disk is presented.     

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.     

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.

     

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

Occurrence of a D-glucosyltransferase in germinating pea seeds and isolation of an endogenous acceptor.

A $\underset{=}{D}$-glucosyltransferase has been found in soluble extracts of germinating pea seeds. The enzyme transfers$\underset{=}{D}$-glucose from uridine diphosphate $\underset{=}{D}$-glucose to a substance which also occurs in the seed extracts. The substance was the only active acceptor of$\underset{=}{D}$-glucose found for this enzymatic reaction. The enzyme has been separated from the acceptor by several purification steps. The acceptor has been purified to near homogeneity. It appears to be a complex glycoside which contains $\underset{=}{L}$-rhamnose, $\underset{=}{D}$-galactose, and $\underset{=}{D}$-glucuronic acid. The aglycone portion of the acceptor has not been characterized, but preliminary findings suggest that it is phenolic in nature.     

A Model of the Quaternary Structure of Enolases, Based on Structural and Evolutionary Analysis of the Octameric Enolase from Bacillus subtilis

Purified enolase from Bacillus subtilis has a native mass of approximately 370 kDa. Since B. subtilis enolase was found to have a subunit mass of 46.58 kDa, the quaternary structure of B. subtilis is octameric. The pl for B. subtilis enolase is 6.1, the pH optimum (pH_o) for activity is 8.1–8.2, and the K _m for 2-PGA is approximately 0.67 mM. Using the dimeric Cα structure of yeast dimeric enolase as a guide, these dimers were arranged as a tetramer of dimers to simulate the electron microscopy image processing obtained for the octameric enolase purified from Thermotoga maritima. This arrangement allowed identification of helix J of one dimer (residues 86–96) and the loop between helix L and strand 1 (HL–S1 loop) of another dimer as possible subunit interaction regions. Alignment of available enolase amino acid sequences revealed that in 16 there are two tandem glycines at the C-terminal end of helix L and the HL–S1 loop is truncated by 4–6 residues relative to the yeast polypeptide, two structural features absent in enolases known to be dimers. From these arrangements and alignments it is proposed that the GG tandem at the C-terminal end of helix L and truncation of the HL–S1 loop may play a critical role in octamer formation of enolases. Interestingly, the sequence features associated with dimeric quaternary structure are found in three phylogenetically disparate groups, suggesting that the ancestral enolase was an octamer and that the dimeric structure has arisen independently multiple times through evolutionary history.     

On the Adjacent Eccentric Distance Sum Index of Graphs

For a given graph G, ε(v) and deg(v) denote the eccentricity and the degree of the vertex v in G, respectively. The adjacent eccentric distance sum index of a graph G is defined as ξsv(G)=∑v∈V(G)ε(v)D(v)deg(v), where D(v)=∑u∈V(G)d(u,v) is the sum of all distances from the vertex v. In this paper we derive some bounds for the adjacent eccentric distance sum index in terms of some graph parameters, such as independence number, covering number, vertex connectivity, chromatic number, diameter and some other graph topological indices.     

pH-, temperature- and ion-dependent oligomerization of Sulfolobus solfataricus recombinant amidase: a study with site-specific mutants

Recombinant amidase from Sulfolobus solfataricus occurred as a dimer of 110 kDa comprising identical subunits. Only dimers were present at pHs above 7.0, but with decreasing pH, dimers associated into octamers, with complete oligomerization occurring at pH 3.0. Oligomerization showed reversible temperature-dependence, with octamer formation increasing with temperature from 36 °C to between 70 and 80° C. Increasing salt concentrations, favored dissociation of the octamers. Among the three investigated factors affecting the dimer–octamer equilibrium, the most important was pH. Among four mutants obtained by site-specific mutagenesis and selection for pH and temperature sensitivity, the T319I and D487N mutant amidases, like that of the native Sulfolobus solfataricus, responded to changes in pH and temperature with a conformational change affecting the dimer–octamer equilibrium. The Y41C and L34P mutant amidases were unaffected by pH and temperature, remaining always in the dimeric state. The differences among mutants in protein conformation must be related to the position of the introduced mutation. Although the L34P and Y41C mutations are located in the helical region 33–48 (LLKLQLESYERLDSLP), which is close to the amino-terminal segment of the protein, the T319I mutation is located in a strand on the surface of the protein, which is far from, and opposite to, the amino-terminal segment. The D487N mutation is located in the center of the protein, far distant from the 33–48 segment. These observations suggest that the segment of the protein closest to the amino-terminus plays a key role in the association of dimers into octamers.     

Rhinos in the Parks: An Island-Wide Survey of the Last Wild Population of the Sumatran Rhinoceros

In the 200 years since the Sumatran rhinoceros was first scientifically described (Fisher 1814), the range of the species has contracted from a broad region in Southeast Asia to three areas on the island of Sumatra and one in Kalimantan, Indonesia. Assessing population and spatial distribution of this very rare species is challenging because of their elusiveness and very low population number. Using an occupancy model with spatial dependency, we assessed the fraction of the total landscape occupied by Sumatran rhinos over a 30,345-km² survey area and the effects of covariates in the areas where they are known to occur. In the Leuser Landscape (surveyed in 2007), the model averaging result of conditional occupancy estimate was ψ^(SE[ψ^])=0.151(0.109) or 2,371.47 km², and the model averaging result of replicated level detection probability p^(SE[p^])=0.252(0.267); in Way Kambas National Park—2008: ψ^(SE[ψ^])=0.468(0.165) or 634.18 km², and p^(SE[p^])=0.138(0.571); and in Bukit Barisan Selatan National Park—2010: ψ^(SE[ψ^])=0.322(0.049) or 819.67 km², and p^(SE[p^])=0.365(0.42). In the Leuser Landscape, rhino occurrence was positively associated with primary dry land forest and rivers, and negatively associated with the presence of a road. In Way Kambas, occurrence was negatively associated with the presence of a road. In Bukit Barisan Selatan, occurrence was negatively associated with presence of primary dryland forest and rivers. Using the probabilities of site occupancy, we developed spatially explicit maps that can be used to outline intensive protection zones for in-situ conservation efforts, and provide a detailed assessment of conserving Sumatran rhinos in the wild. We summarize our core recommendation in four points: consolidate small population, strong protection, determine the percentage of breeding females, and recognize the cost of doing nothing. To reduce the probability of poaching, here we present only the randomized location of site level occupancy in our result while retaining the overall estimation of occupancy for a given area.