Computationally restoring the potency of a clinical antibody against Omicron.

The COVID-19 pandemic underscored the promise of monoclonal antibody-based prophylactic and therapeutic drugs1-3 and revealed how quickly viral escape can curtail effective options4,5. When the SARS-CoV-2 Omicron variant emerged in 2021, many antibody drug products lost potency, including Evusheld and its constituent, cilgavimab4-6. Cilgavimab, like its progenitor COV2-2130, is a class 3 antibody that is compatible with other antibodies in combination4 and is challenging to replace with existing approaches. Rapidly modifying such high-value antibodies to restore efficacy against emerging variants is a compelling mitigation strategy. We sought to redesign and renew the efficacy of COV2-2130 against Omicron BA.1 and BA.1.1 strains while maintaining efficacy against the dominant Delta variant. Here we show that our computationally redesigned antibody, 2130-1-0114-112, achieves this objective, simultaneously increases neutralization potency against Delta and subsequent variants of concern, and provides protection in vivo against the strains tested: WA1/2020, BA.1.1 and BA.5. Deep mutational scanning of tens of thousands of pseudovirus variants reveals that 2130-1-0114-112 improves broad potency without increasing escape liabilities. Our results suggest that computational approaches can optimize an antibody to target multiple escape variants, while simultaneously enriching potency. Our computational approach does not require experimental iterations or pre-existing binding data, thus enabling rapid response strategies to address escape variants or lessen escape vulnerabilities.

Three-dimensional high-definition exoscope (Kestrel View II) in anterior cervical discectomy and fusion: a valid alternative to operative microscope-assisted surgery.

BACKGROUND: Operative microscope (OM) has greatly advanced modern spine surgery, but remains limited by several drawbacks. Therefore, a three-dimensional (3D) high-definition (HD) exoscope (EX) (Kestrel View II, Mataka Kohli, Japan) system has been developed and used as an alternative to the OM. The aim of this study was to assess and compare the perioperative data and clinical outcomes of anterior cervical discectomy and fusion (ACDF) procedure with either an EX or OM. METHODS: Forty-eight patients with cervical spondylotic myelopathy (CSM) underwent ACDF assisted by the EX or OM between January 2019 and December 2019. We collected and compared data on operative time, intraoperative bleeding, postoperative hospitalization stay, complications, and clinical outcomes between the two groups. The clinical outcomes were evaluated by using visual analogue scale (VAS) scores, Japanese Orthopedic Association (JOA) scores, the recovery rate of JOA scores, and Odom criteria. RESULTS: The operative time in the EX group was significantly shorter than that in the OM group (P < 0.05). The VAS and JOA scores were significantly improved in both groups after surgery (P < 0.05). In addition, the VAS scores in the EX group were significantly lower than those in the OM group at 1 week postoperatively (P < 0.05). The good-to-excellent outcome rates were 90.48 and 88.89% in the EX group and OM group, respectively, whereas the complication occurrence rates of the EX group and OM group were 4.76 and 11.11%, respectively. CONCLUSIONS: EX-assisted and OM-assisted ACDF resulted in similar clinical outcomes for CSM, while EX-assisted surgery may be related to a short operative time and fewer complications.

Binding Affinity Prediction by Pairwise Function Based on Neural Network.

We present a new approach to estimate the binding affinity from given three-dimensional poses of protein-ligand complexes. In this scheme, every protein-ligand atom pair makes an additive free-energy contribution. The sum of these pairwise contributions then gives the total binding free energy or the logarithm of the dissociation constant. The pairwise contribution is calculated by a function implemented via a neural network that takes the properties of the two atoms and their distance as input. The pairwise function is trained using a portion of the PDBbind 2018 data set. The model achieves good accuracy for affinity predictions when evaluated with PDBbind 2018 and with the CASF-2016 benchmark, comparing favorably to many scoring functions such as that of AutoDock Vina. The framework here may be extended to incorporate other factors to further improve its accuracy and power.

Calculating Single-Channel Permeability and Conductance from Transition Paths.

Permeability and conductance are the major transport properties of membrane channels, quantifying the rate of channel crossing by the solute. It is highly desirable to calculate these quantities in all-atom molecular dynamics simulations. When the solute crossing rate is low, however, direct methods would require prohibitively long simulations, and one thus typically adopts alternative strategies based on the free energy of single solute along the channel. Here we present a new method to calculate the crossing rate by initiating unbiased trajectories in which the solute is released at the free energy barrier. In this method, the total time the solute spends in the barrier region during a channel crossing (transition path) is used to determine the kinetic rate. Our method achieves a significantly higher statistical accuracy than the classical reactive flux method, especially for diffusive barrier crossing. Our test on ion permeation through a carbon nanotube verifies that the method correctly predicts the crossing rate and reproduces the spontaneous crossing events as in long equilibrium simulations. The rigorous and efficient method here will be valuable for quantitatively connecting simulations to experimental measurement of membrane channels.

Structural Model of the Tubular Assembly of the Rous Sarcoma Virus Capsid Protein.

The orthoretroviral capsid protein (CA) assembles into polymorphic capsids, whose architecture, assembly, and stability are still being investigated. The N-terminal and C-terminal domains of CA (NTD and CTD, respectively) engage in both homotypic and heterotypic interactions to create the capsid. Hexameric turrets formed by the NTD decorate the majority of the capsid surface. We report nearly complete solid-state NMR (ssNMR) resonance assignments of Rous sarcoma virus (RSV) CA, assembled into hexamer tubes that mimic the authentic capsid. The ssNMR assignments show that, upon assembly, large conformational changes occur in loops connecting helices, as well as the short 310 helix initiating the CTD. The interdomain linker becomes statically disordered. Combining constraints from ssNMR and cryo-electron microscopy (cryo-EM), we establish an atomic resolution model of the RSV CA tubular assembly using molecular dynamics flexible fitting (MDFF) simulations. On the basis of comparison of this MDFF model with an earlier-derived crystallographic model for the planar assembly, the induction of curvature into the RSV CA hexamer lattice arises predominantly from reconfiguration of the NTD-CTD and CTD trimer interfaces. The CTD dimer and CTD trimer interfaces are also intrinsically variable. Hence, deformation of the CA hexamer lattice results from the variable displacement of the CTDs that surround each hexameric turret. Pervasive H-bonding is found at all interdomain interfaces, which may contribute to their malleability. Finally, we find helices at the interfaces of HIV and RSV CA assemblies have very different contact angles, which may reflect differences in the capsid assembly pathway for these viruses.

Parameter Optimization for Interaction between C-Terminal Domains of HIV-1 Capsid Protein.

HIV-1 capsid proteins (CAs) assemble into a capsid that encloses the viral RNA. The binding between a pair of C-terminal domains (CTDs) constitutes a major interface in both the CA dimers and the large CA assemblies. Here, we attempt to use a general residue-level coarse-grained model to describe the interaction between two isolated CTDs in Monte Carlo simulations. With the standard parameters that depend only on the residue types, the model predicts a much weaker binding in comparison to the experiments. Detailed analysis reveals that some Lennard-Jones parameters are not compatible with the experimental CTD dimer structure, thus resulting in an unfavorable interaction energy. To improve the model for the CTD binding, we introduce ad hoc modifications to a small number of Lennard-Jones parameters for some specific pairs of residues at the binding interface. Through a series of extensive Monte Carlo simulations, we identify the optimal parameters for the CTD-CTD interactions. With the refined model parameters, both the binding affinity (with a dissociation constant of 13 ± 2 µM) and the binding mode are in good agreement with the experimental data. This study demonstrates that the general interaction model based on the Lennard-Jones potential, with some modest adjustment of the parameters for key residues, could correctly reproduce the reversible protein binding, thus potentially applicable for simulating the thermodynamics of the CA assemblies.

Curcumin Suppresses Epithelial-Mesenchymal Transition of Renal Tubular Epithelial Cells through the Inhibition of Akt/mTOR Pathway.

Curcumin has exhibited a protective effect against development of renal fibrosis in animal models, however, its underlying molecular mechanisms are largely unclear. Therefore, we investigated the anti-fibrosis effects of curcumin in transforming growth factor-ß1 (TGF-ß1)-induced epithelial-to-mesenchymal transition (EMT), and the mechanism by which it mediates the phosphatidylinositol 3-kinase (PI3K)/Akt pathway. Human kidney tubular epithelial cells (HKCs) were treated with TGF-ß1 or curcumin alone, or TGF-ß1 in combination with curcumin. The effect of curcumin on cell proliferation was measured by 3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide (MTT) assay. Expression of E-cadherin, cytokeratin, vimentin, alpha smooth muscle actin (α-SMA), fibroblast-specific protein 1 (FSP1) and key proteins of Akt/mammalian target of rapamycin (mTOR) pathway were analyzed by immunocytochemistry, real-time PCR and Western blot. Low dose curcumin (3.125 and 25 µmol/L) effectively promoted HKC proliferation. When HKCs were co-incubated with TGF-ß1 and curcumin for 72 h, curcumin maintained the epithelial morphology in a dose-dependent manner, decreased expression of vimentin, α-SMA and FSP1 normally induced by TGF-ß1, and increased expression of E-cadherin, cytokeratin. Importantly, we found that curcumin reduced Akt, mTOR and P70S6K phosphorylation, effectively suppressing the activity of the Akt/mTOR pathway in HKCs. Curcumin also promoted HKC proliferation, and antagonized TGF-ß1-driven EMT through the inhibition of Akt/mTOR pathway activity, which may suggest an alternative therapy for renal fibrosis.

Calculating transition rates from durations of transition paths.

Calculating the kinetic rates for rare transitions is an important objective for molecular simulations. Here I prove equalities that relate the transition rates to the equilibrium free energy and the statistics of the transition paths. In particular, the durations of the transition paths within given intervals of the reaction coordinate provide the kinetic pre-factor in the rate formula. Based on the available free energy, the transition rates can further be rigorously calculated by initiating forward and backward simulations and evaluating the duration of each transition path. Validation on a model system confirms that the approach correctly predicts the transition rates from the simulations and demonstrates that whereas the relations here are general and valid for any chosen reaction coordinate, a good reaction coordinate will enable a more efficient sampling of the transition paths and thus a more reliable rate calculation.

Mechanism of polymorphism and curvature of HIV capsid assemblies probed by 3D simulations with a novel coarse grain model.

BACKGROUND: During the maturation process, HIV capsid proteins self-assemble into polymorphic capsids. The strong polymorphism precludes high resolution structural characterization under in vivo conditions. In spite of the determination of structural models for various in vitro assemblies of HIV capsid proteins, the assembly mechanism is still not well-understood. METHODS: We report 3D simulations of HIV capsid proteins by a novel coarse grain model that captures the backbone of the rigid segments in the protein accurately. The effects of protein dynamics on assembly are emulated by a static ensemble of subunits in conformations derived from molecular dynamics simulation. RESULTS: We show that HIV capsid proteins robustly assemble into hexameric lattices in a range of conditions where trimers of dimeric subunits are the dominant oligomeric intermediates. Variations of hexameric lattice curvatures are observed in simulations with subunits of variable inter-domain orientations mimicking the conformation distribution in solution. Simulations with subunits based on pentameric structural models lead to assembly of sharp curved structures resembling the tips of authentic HIV capsids, along a distinct pathway populated by tetramers and pentamers with the characteristic quasi-equivalency of viral capsids. CONCLUSIONS: Our results suggest that the polymorphism assembly is triggered by the inter-domain dynamics of HIV capsid proteins in solution. The assembly of highly curved structures arises from proteins in conformation with a highly specific inter-domain orientation. SIGNIFICANCE: Our work proposes a mechanism of HIV capsid assembly based on available structural data, which can be readily verified. Our model can be applied to other large biomolecular assemblies.

Monte Carlo Simulations of HIV Capsid Protein Homodimer.

Capsid protein (CA) is the building block of virus coats. To help understand how the HIV CA proteins self-organize into large assemblies of various shapes, we aim to computationally evaluate the binding affinity and interfaces in a CA homodimer. We model the N- and C-terminal domains (NTD and CTD) of the CA as rigid bodies and treat the five-residue loop between the two domains as a flexible linker. We adopt a transferrable residue-level coarse-grained energy function to describe the interactions between the protein domains. In seven extensive Monte Carlo simulations with different volumes, a large number of binding/unbinding transitions between the two CA proteins are observed, thus allowing a reliable estimation of the equilibrium probabilities for the dimeric vs monomeric forms. The obtained dissociation constant for the CA homodimer from our simulations, 20-25 µM, is in reasonable agreement with experimental measurement. A wide range of binding interfaces, primarily between the NTDs, are identified in the simulations. Although some observed bound structures here closely resemble the major binding interfaces in the capsid assembly, they are statistically insignificant in our simulation trajectories. Our results suggest that although the general purpose energy functions adopted here could reasonably reproduce the overall binding affinity for the CA homodimer, further adjustment would be needed to accurately represent the relative strength of individual binding interfaces.

Role of electrostatic interactions in binding of peptides and intrinsically disordered proteins to their folded targets. 1. NMR and MD characterization of the complex between the c-Crk N-SH3 domain and the peptide Sos.

Intrinsically disordered proteins (IDPs) often rely on electrostatic interactions to bind their structured targets. To obtain insight into the mechanism of formation of the electrostatic encounter complex, we investigated the binding of the peptide Sos (PPPVPPRRRR), which serves as a minimal model for an IDP, to the c-Crk N-terminal SH3 domain. Initially, we measured ¹5N relaxation rates at two magnetic field strengths and determined the binding shifts for the complex of Sos with wild-type SH3. We have also recorded a 3 µs molecular dynamics (MD) trajectory of this complex using the Amber ff99SB*-ILDN force field. The comparison of the experimental and simulated data shows that MD simulation consistently overestimates the strength of salt bridge interactions at the binding interface. The series of simulations using other advanced force fields also failed to produce any satisfactory results. To address this issue, we have devised an empirical correction to the Amber ff99SB*-ILDN force field whereby the Lennard-Jones equilibrium distance for the nitrogen-oxygen pair across the Arg-to-Asp and Arg-to-Glu salt bridges has been increased by 3%. Implementing this correction resulted in a good agreement between the simulations and the experiment. Adjusting the strength of salt bridge interactions removed a certain amount of strain contained in the original MD model, thus improving the binding of the hydrophobic N-terminal portion of the peptide. The arginine-rich C-terminal portion of the peptide, freed from the effect of the overstabilized salt bridges, was found to interconvert more rapidly between its multiple conformational states. The modified MD protocol has also been successfully used to simulate the entire binding process. In doing so, the peptide was initially placed high above the protein surface. It then arrived at the correct bound pose within ∼2 Å of the crystallographic coordinates. This simulation allowed us to analyze the details of the dynamic binding intermediate, i.e., the electrostatic encounter complex. However, an experimental characterization of this transient, weakly populated state remains out of reach. To overcome this problem, we designed the double mutant of c-Crk N-SH3 in which mutations Y186L and W169F abrogate tight Sos binding and shift the equilibrium toward the intermediate state resembling the electrostatic encounter complex. The results of the combined NMR and MD study of this engineered system will be reported in the next part of this paper.

Hexagonal Lattices of HIV Capsid Proteins Explored by Simulations Based on a Thermodynamically Consistent Model.

HIV capsid proteins (CAs) may self-assemble into a variety of shapes under in vivo and in vitro conditions. Here, we employed simulations based on a residue-level coarse-grained (CG) model with full conformational flexibility to investigate hexagonal lattices, which are the underlying structural pattern for CA aggregations. Facilitated by enhanced sampling simulations to rigorously calculate CA dimerization and polymerization affinities, we calibrated our model to reproduce the experimentally measured affinities. Using the calibrated model, we performed unbiased simulations on several large systems consisting of 1512 CA subunits, allowing reversible binding and unbinding of the CAs in a thermodynamically consistent manner. In one simulation, a preassembled hexagonal CA sheet developed spontaneous curvatures reminiscent of those observed in experiments, and the edges of the sheet exhibited local curvatures larger than those of the interior. In other simulations starting with randomly distributed CAs at different concentrations, existing CA assemblies grew by binding free capsomeres to the edges and by merging with other assemblies. At high CA concentrations, rapid establishment of predominant aggregates was followed by much slower adjustments toward more regular hexagonal lattices, with increasing numbers of intact CA hexamers and pentamers being formed. Our approach of adapting a general CG model to specific systems by using experimental binding data represents a practical and effective strategy for simulating and elucidating intricate protein aggregations.

Collective diffusion model for ion conduction through microscopic channels.

Ion conduction through microscopic channels is of central importance in both biology and nanotechnology. To better understand the current-voltage (I-V) dependence of ion channels, here we describe and prove a collective diffusion model that quantitatively relates the spontaneous ion permeation at equilibrium to the stationary ionic fluxes driven by small voltages. The model makes it possible to determine the channel conductance in the linear I-V range from equilibrium simulations without the application of a voltage. To validate the theory, we perform molecular-dynamics simulations on two channels-a conical-shaped nanopore and the transmembrane pore of an α-hemolysin-under both equilibrium and nonequilibrium conditions. The simulations reveal substantial couplings between the motions of cations and anions, which are effectively captured by the collective coordinate in the model. Although the two channels exhibit very different linear ranges in the I-V curves, in both cases the channel conductance at small voltages is in reasonable agreement with the prediction from the equilibrium simulation. The simulations also suggest that channel charges, rather than geometric asymmetry, play a more prominent role in current rectification.

Stromal cell-derived factor-1 enhances wound healing through recruiting bone marrow-derived mesenchymal stem cells to the wound area and promoting neovascularization.

Stromal cell-derived factor-1 (SDF-1) is a potent chemokine for bone marrow-derived stromal stem cells (BMSCs) that express CXCR4, the receptor for SDF-1. SDF-1 is considered to play an important role in the trafficking of BMSCs. We investigated the contribution of SDF-1 to the recruitment of BMSCs to the wound area and its promotion of wound repair and neovascularization. BMSCs were pretreated with or without anti-CXCR4 blocking antibody and combined with CM-DiI label, and injected via the tail vein into mice with full-thickness skin wounds on the dorsum. Simultaneously, anti-SDF-1 antibody was injected into local wounds in another group of mice. The results show that blockade of CXCR4 on either infused BMSCs or SDF-1 in the host wounds (1) dramatically impaired the number of infused BMSCs being recruited to the injured tissue, (2) reduced the expression of growth factors involved in the repair of injured tissue such as vascular endothelial growth factor, basic fibroblast growth factor and transforming growth factor beta 1, (3) decreased the resultant neovascularization, and (4) retarded wound healing. Taken together, the findings indicate that the SDF-1/CXCR4 signal pathway facilitates wound healing through augmenting BMSC recruitment to wound tissues, responsive secretion of growth factors by BMSCs and neovascularization in the wound area.

Pore opening and closing of a pentameric ligand-gated ion channel.

Nerve signaling in humans and chemical sensing in bacteria both rely on the controlled opening and closing of the ion-conducting pore in pentameric ligand-gated ion channels. With the help of a multiscale simulation approach that combines mixed elastic network model calculations with molecular dynamics simulations, we study the opening and closing of the pore in Gloeobacter violaceus channel GLIC at atomic resolution. In our simulations of the GLIC transmembrane domain, we first verify that the two endpoints of the transition are open and closed to sodium ion conduction, respectively. We then show that a two-stage tilting of the pore-lining helices induces cooperative drying and iris-like closing of the channel pore. From the free energy profile of the gating transition and from unrestrained simulations, we conclude that the pore of the isolated GLIC transmembrane domain closes spontaneously. The mechanical work of opening the pore is performed primarily on the M2-M3 loop. Strong interactions of this short and conserved loop with the extracellular domain are therefore crucial to couple ligand binding to channel opening.

Three-dimensional High-definition Exoscope in Minimally Invasive Transforaminal Lumbar Interbody Fusion: A Retrospective Cohort Study.

OBJECTIVES: The operative microscope (OM) has revolutionized the field of modern spine surgery, however, it remains limited by several drawbacks. Recently, the exoscope (EX) system has been designed to assistant spine surgery. It provides a three-dimensional (3D) high-definition (HD) operative experience and becomes an alternative to the OM. The aim of the study was to evaluate the clinical outcomes, advantages and limitations of EX-assisted minimally invasive transforaminal lumbar interbody fusion (EMIS-TLIF) and OM-assisted MIS-TLIF (OMIS-TLIF). METHODS: The clinical outcomes were assessed in 47 patients with lumbar degenerative diseases (LDD) who underwent MIS-TLIF assisted with the OM or EX between January 2019 and September 2020. A total of 22 were treated with EMIS-TLIF, and 25 received OMIS-TLIF. Perioperative parameters (including sex, age, number of fusion levels and body mass index), perioperative parameters (operation time, intraoperative blood loss, postoperative drainage, postoperative hospitalization stay, and duration of follow-up), visual analogue scale (VAS) of back pain, VAS of leg pain, Oswestry disability index (ODI) scores and clinical outcomes were assessed and compared. Image quality, handling of equipment, ergonomics, 3D glasses and educational usefulness were scored according to a questionnaire. RESULTS: Operation time in the OMIS-TLIF group (121.92 ± 16.92 min) was significantly increased compared with that in the EMIS-TLIF group (111.00 ± 19.87 min) (P < 0.05). The VAS of the back pain and ODI scores in the EMIS-TLIF group were significantly lower compared with the OMIS-TLIF group at 1 week postoperatively (P < 0.05). The good-excellent outcomes rate was 90.91% in the EMIS-TLIF group and 88.00% in the OMIS-TLIF group, and there was no significant difference. A total of 44 visits completed the questionnaire. The results of the questionnaire showed that the EX has exhibited advantages regarding handing of equipment, ergonomics and educational usefulness, and comparable image quality as compared with the OM, however, operating surgeons complained uncomfortable sensation when wearing 3D glasses. CONCLUSIONS: The EMIS-TLIF was a safe and effective procedure in the management of LDD as compared with the OMIS-LIF. Meanwhile, EMIS-TLIF might resulted in a short operation time.

Computationally restoring the potency of a clinical antibody against SARS-CoV-2 Omicron subvariants.

The COVID-19 pandemic underscored the promise of monoclonal antibody-based prophylactic and therapeutic drugs1-3, but also revealed how quickly viral escape can curtail effective options4,5. With the emergence of the SARS-CoV-2 Omicron variant in late 2021, many clinically used antibody drug products lost potency, including Evusheld™ and its constituent, cilgavimab4,6. Cilgavimab, like its progenitor COV2-2130, is a class 3 antibody that is compatible with other antibodies in combination4 and is challenging to replace with existing approaches. Rapidly modifying such high-value antibodies with a known clinical profile to restore efficacy against emerging variants is a compelling mitigation strategy. We sought to redesign COV2-2130 to rescue in vivo efficacy against Omicron BA.1 and BA.1.1 strains while maintaining efficacy against the contemporaneously dominant Delta variant. Here we show that our computationally redesigned antibody, 2130-1-0114-112, achieves this objective, simultaneously increases neutralization potency against Delta and many variants of concern that subsequently emerged, and provides protection in vivo against the strains tested, WA1/2020, BA.1.1, and BA.5. Deep mutational scanning of tens of thousands pseudovirus variants reveals 2130-1-0114-112 improves broad potency without incurring additional escape liabilities. Our results suggest that computational approaches can optimize an antibody to target multiple escape variants, while simultaneously enriching potency. Because our approach is computationally driven, not requiring experimental iterations or pre-existing binding data, it could enable rapid response strategies to address escape variants or pre-emptively mitigate escape vulnerabilities.

Drying transition in the hydrophobic gate of the GLIC channel blocks ion conduction.

The theoretical prediction of water drying transitions near nonpolar surfaces has stimulated an intensive search for biological processes exploiting this extreme form of hydrophobicity. Here we quantitatively demonstrate that drying of a hydrophobic constriction is the major determinant of ion conductance in the GLIC pentameric ion channel. Molecular-dynamics simulations show that in the closed state, the channel conductance is ∼12 orders-of-magnitude lower than in the open state. This large drop in conductance is remarkable because even in the functionally closed conformation the pore constriction remains wide enough for the passage of sodium ions, aided by a continuous bridge of ∼12 water molecules. However, we find that the free energy cost of hydrating the hydrophobic gate is large, accounting almost entirely for the energetic barrier blocking ion passage. The free energies of transferring a sodium ion into a prehydrated gate in functionally closed and open states differ by only 1.2 kcal/mol, compared to an 11 kcal/mol difference in the costs of hydrating the hydrophobic gate. Conversely, ion desolvation effects play only minor roles in GLIC ion channel gating. Our simulations help rationalize experiments probing the gating kinetics of the nicotinic acetylcholine receptor in response to mutations of pore-lining residues. The molecular character and phase behavior of water should thus be included in quantitative descriptions of ion channel gating.

Convergence and error estimation in free energy calculations using the weighted histogram analysis method.

The weighted histogram analysis method (WHAM) has become the standard technique for the analysis of umbrella sampling simulations. In this article, we address the challenges (1) of obtaining fast and accurate solutions of the coupled nonlinear WHAM equations, (2) of quantifying the statistical errors of the resulting free energies, (3) of diagnosing possible systematic errors, and (4) of optimally allocating of the computational resources. Traditionally, the WHAM equations are solved by a fixed-point direct iteration method, despite poor convergence and possible numerical inaccuracies in the solutions. Here, we instead solve the mathematically equivalent problem of maximizing a target likelihood function, by using superlinear numerical optimization algorithms with a significantly faster convergence rate. To estimate the statistical errors in one-dimensional free energy profiles obtained from WHAM, we note that for densely spaced umbrella windows with harmonic biasing potentials, the WHAM free energy profile can be approximated by a coarse-grained free energy obtained by integrating the mean restraining forces. The statistical errors of the coarse-grained free energies can be estimated straightforwardly and then used for the WHAM results. A generalization to multidimensional WHAM is described. We also propose two simple statistical criteria to test the consistency between the histograms of adjacent umbrella windows, which help identify inadequate sampling and hysteresis in the degrees of freedom orthogonal to the reaction coordinate. Together, the estimates of the statistical errors and the diagnostics of inconsistencies in the potentials of mean force provide a basis for the efficient allocation of computational resources in free energy simulations.

Large-scale application of free energy perturbation calculations for antibody design.

Alchemical free energy perturbation (FEP) is a rigorous and powerful technique to calculate the free energy difference between distinct chemical systems. Here we report our implementation of automated large-scale FEP calculations, using the Amber software package, to facilitate antibody design and evaluation. In combination with Hamiltonian replica exchange, our FEP simulations aim to predict the effect of mutations on both the binding affinity and the structural stability. Importantly, we incorporate multiple strategies to faithfully estimate the statistical uncertainties in the FEP results. As a case study, we apply our protocols to systematically evaluate variants of the m396 antibody for their conformational stability and their binding affinity to the spike proteins of SARS-CoV-1 and SARS-CoV-2. By properly adjusting relevant parameters, the particle collapse problems in the FEP simulations are avoided. Furthermore, large statistical errors in a small fraction of the FEP calculations are effectively reduced by extending the sampling, such that acceptable statistical uncertainties are achieved for the vast majority of the cases with a modest total computational cost. Finally, our predicted conformational stability for the m396 variants is qualitatively consistent with the experimentally measured melting temperatures. Our work thus demonstrates the applicability of FEP in computational antibody design.