Machine Learning Applied to Determine the Molecular Descriptors Responsible for the Viscosity Behavior of Concentrated Therapeutic Antibodies.

Predicting the solution viscosity of monoclonal antibody (mAb) drug products remains as one of the main challenges in antibody drug design, manufacturing, and delivery. In this work, the concentration-dependent solution viscosity of 27 FDA-approved mAbs was measured at pH 6.0 in 10 mM histidine-HCl. Six mAbs exhibited high viscosity (>30 cP) in solutions at 150 mg/mL mAb concentration. Combining molecular modeling and machine learning feature selection, we found that the net charge in the mAbs and the amino acid composition in the Fv region are key features which govern the viscosity behavior. For mAbs whose behavior was not dominated by charge effects, we observed that high viscosity is correlated with more hydrophilic and fewer hydrophobic residues in the Fv region. A predictive model based on the net charges of mAbs and a high viscosity index is presented as a fast screening tool for classifying low- and high-viscosity mAbs.

Identification and elucidation of proline-rich antimicrobial peptides with enhanced potency and delivery.

Proline-rich antimicrobial peptides (PrAMPs) kill bacteria via a nonlytic mechanism in which they permeate through the outer membrane, utilize protein-mediated transport across the inner membrane, and target the ribosome to inhibit protein synthesis. We previously reported that substitutions of oncocin ( VDKPPYLPRPRPPRRIYNR-NH2 ) with a pair of cationic residues improved the antimicrobial activity. In this study, we applied the design protocol to three other PrAMPs: apidaecin-1b, pyrrhocoricin, and bactenecin 7(1-16) and found that the substitutions (R4K and I8K/R) for apidaecin-1b improve the activity by twofold (p < .05) against nonpathogenic Escherichia coli. Moreover, the substitutions (L7K/R and R14K) for pyrrhocoricin improve the activity by 2-10-fold (p < .05) against some strains of E. coli and Salmonella Typhimurium. We also performed activity tests against inner membrane protein (SbmA or YgdD) knockout strains. The result is consistent with previous studies that SbmA is the major transporter for apidaecin-1b and pyrrhocoricin derivatives. However, bactenecin 7(1-16) functions independently of these transporters. In addition, several apidaecin-1b derivatives exhibit enhanced activity relative to wild-type only in the absence of SbmA, which is consistent with mutations that enhance transport across the inner membrane. A high performance liquid chromatography-based kinetic assay for cellular association and internalization demonstrates that the selected cationic mutations can improve cellular association in minimal media, but this enhanced association is not required for increased activity, which suggests the importance of inner membrane transport. These functional studies on cationic mutants of PrAMPs advance understanding of potency and mechanism and advance the ability to engineer improved antimicrobials as evidenced by the identification of the pyrrhocoricin mutant (L7R and R14K) with 10-fold elevated potency against pathogenic E. coli.

Internal coordinate density of state from molecular dynamics simulation.

The vibrational density of states (DoS), calculated from the Fourier transform of the velocity autocorrelation function, provides profound information regarding the structure and dynamic behavior of a system. However, it is often difficult to identify the exact vibrational mode associated with a specific frequency if the DoS is determined based on velocities in Cartesian coordinates. Here, the DoS is determined based on velocities in internal coordinates, calculated from Cartesian atomic velocities using a generalized Wilson's B-matrix. The DoS in internal coordinates allows for the correct detection of free dihedral rotations that may be mistaken as hindered rotation in Cartesian DoS. Furthermore, the pronounced enhancement of low frequency modes in Cartesian DoS for macromolecules should be attributed to the coupling of dihedral and angle motions. The internal DoS, thus deconvolutes the internal motions and provides fruitful insights to the dynamic behaviors of a system.

DeepSP: Deep learning-based spatial properties to predict monoclonal antibody stability.

Therapeutic antibody development faces challenges due to high viscosities and aggregation tendencies. The spatial charge map (SCM) and spatial aggregation propensity (SAP) are computational techniques that aid in predicting viscosity and aggregation, respectively. These methods rely on structural data derived from molecular dynamics (MD) simulations, which are computationally demanding. DeepSCM, a deep learning surrogate model based on sequence information to predict SCM, was recently developed to screen high-concentration antibody viscosity. This study further utilized a dataset of 20,530 antibody sequences to train a convolutional neural network deep learning surrogate model called Deep Spatial Properties (DeepSP). DeepSP directly predicts SAP and SCM scores in different domains of antibody variable regions based solely on their sequences without performing MD simulations. The linear correlation coefficient between DeepSP scores and MD-derived scores for 30 properties achieved values between 0.76 and 0.96 with an average of 0.87. DeepSP descriptors were employed as features to build machine learning models to predict the aggregation rate of 21 antibodies, and the performance is similar to the results obtained from the previous study using MD simulations. This result demonstrates that the DeepSP approach significantly reduces the computational time required compared to MD simulations. The DeepSP model enables the rapid generation of 30 structural properties that can also be used as features in other research to train machine learning models for predicting various antibody stability using sequences only. DeepSP is freely available as an online tool via https://deepspwebapp.onrender.com and the codes and parameters are freely available at https://github.com/Lailabcode/DeepSP.

Teaching biologics design using molecular modeling and simulations.

Teaching chemistry and biology students about biologics design remains challenging despite its increasing importance in pharmaceutical development. Monoclonal antibodies, commonly called mAbs, are the most popular biologics. They have been developed into drugs to treat various diseases in the past decades. Multiple challenges exist for designing proper formulations to stabilize mAbs, such as preventing aggregation and mitigating viscosity. Molecular modeling and simulations can improve pharmaceutical products by examining the interactions between mAbs and other compounds, such as excipients. To introduce students to biopharmaceuticals, eight students at the Stevens Institute of Technology participated in a semester-long course to learn the challenges of pharmaceutical development and different computational skills to study biologics design. The students started with a limited background in this field. Throughout one semester, they were introduced to various literature and software tools for modeling antibodies and studying their interactions with excipients. This paper aims to develop a course structure to be replicated at other universities and institutions to teach biopharmaceutical development to students.

Determination of Metabolic Fluxes by Deep Learning of Isotope Labeling Patterns.

Fluxomics offers a direct readout of metabolic state but relies on indirect measurement. Stable isotope tracers imprint flux-dependent isotope labeling patterns on metabolites we measure; however, the relationship between labeling patterns and fluxes remains elusive. Here we innovate a two-stage machine learning framework termed ML-Flux that streamlines metabolic flux quantitation from isotope tracing. We train machine learning models by simulating atom transitions across five universal metabolic models starting from 26 13C-glucose, 2H-glucose, and 13C-glutamine tracers within feasible flux space. ML-Flux employs deep-learning-based imputation to take variable measurements of labeling patterns as input and successive neural networks to convert the ensuing comprehensive labeling information into metabolic fluxes. Using ML-Flux with multi-isotope tracing, we obtain fluxes through central carbon metabolism that are comparable to those from a least-squares method but orders-of-magnitude faster. ML-Flux is deployed as a webtool to expand the accessibility of metabolic flux quantitation and afford actionable information on metabolism.

Rapid determination of entropy and free energy of mixtures from molecular dynamics simulations with the two-phase thermodynamic model.

The two-phase thermodynamic (2PT) model is generalized to determine the thermodynamic properties of mixtures. In this method, the vibrational density of states (DoS), obtained from the Fourier transform of the velocity autocorrelation function, and quantum statistics are combined to determine the entropy and free energy from the trajectory of a molecular dynamics simulation. In particular, the calculated DoS is decomposed into a solid-like and a gas-like component through the fluidicity parameter, allowing for treatments for the anharmonic effects in fluids. The 2PT method has been shown to provide reliable thermodynamic properties of pure substances over the whole phase diagram with only about a 20 ps MD trajectory. Here we show how the 2PT method can be used for mixtures with the same degree of accuracy and efficiency. We have examined the 2PT determined excess Gibbs free energies of Lennard-Jones (LJ) mixtures over a wide range of conditions (1 ≤ T* ≤ 3, 0.5 ≤ P* ≤ 2.5, 1 ≤ σ(BB)/σ(AA) ≤ 2, and 1 ≤ ε(BB)/ε(AA) ≤ 2), including the change of the off-diagonal LJ interactions. The 2PT determined values are in good agreement with those from Widom insertion or thermodynamic integration (TI). Our results suggest that the 2PT method can be a powerful method for understanding thermodynamic properties in more complicated multicomponent systems.

DeepSCM: An efficient convolutional neural network surrogate model for the screening of therapeutic antibody viscosity.

Predicting high concentration antibody viscosity is essential for developing subcutaneous administration. Computer simulations provide promising tools to reach this aim. One such model is the spatial charge map (SCM) proposed by Agrawal and coworkers (mAbs. 2015, 8(1):43-48). SCM applies molecular dynamics simulations to calculate a score for the screening of antibody viscosity at high concentrations. However, molecular dynamics simulations are computationally costly and require structural information, a significant application bottleneck. In this work, high throughput computing was performed to calculate the SCM scores for 6596 nonredundant antibody variable regions. A convolutional neural network surrogate model, DeepSCM, requiring only sequence information, was then developed based on this dataset. The linear correlation coefficient of the DeepSCM and SCM scores achieved 0.9 on the test set (N = 1320). The DeepSCM model was applied to screen the viscosity of 38 therapeutic antibodies that SCM correctly classified and resulted in only one misclassification. The DeepSCM model will facilitate high concentration antibody viscosity screening. The code and parameters are freely available at https://github.com/Lailabcode/DeepSCM.

Machine learning prediction of antibody aggregation and viscosity for high concentration formulation development of protein therapeutics.

Machine learning has been recently used to predict therapeutic antibody aggregation rates and viscosity at high concentrations (150 mg/ml). These works focused on commercially available antibodies, which may have been optimized for stability. In this study, we measured accelerated aggregation rates at 45°C and viscosity at 150 mg/ml for 20 preclinical and clinical-stage antibodies. Features obtained from molecular dynamics simulations of the full-length antibody and sequences were used for machine learning model construction. We found a k-nearest neighbors regression model with two features, spatial positive charge map on the CDRH2 and solvent-accessible surface area of hydrophobic residues on the variable fragment, gives the best performance for predicting antibody aggregation rates (r = 0.89). For the viscosity classification model, the model with the highest accuracy is a logistic regression model with two features, spatial negative charge map on the heavy chain variable region and spatial negative charge map on the light chain variable region. The accuracy and the area under precision recall curve of the classification model from validation tests are 0.86 and 0.70, respectively. In addition, we combined data from another 27 commercial mAbs to develop a viscosity predictive model. The best model is a logistic regression model with two features, number of hydrophobic residues on the light chain variable region and net charges on the light chain variable region. The accuracy and the area under precision recall curve of the classification model are 0.85 and 0.6, respectively. The aggregation rates and viscosity models can be used to predict antibody stability to facilitate pharmaceutical development.

Calculation of therapeutic antibody viscosity with coarse-grained models, hydrodynamic calculations and machine learning-based parameters.

High viscosity presents a challenge for manufacturing and drug delivery of therapeutic antibodies. The viscosity is determined by protein-protein interactions among many antibodies. Molecular simulation is a promising method to study protein-protein interactions; however, all-atom models do not allow the simulation of multiple molecules, which is necessary to compute viscosity directly. Coarse-grained models, on the other hand can do this. In this work, a 12-bead coarse-grained model based on Swan and coworkers (J. Phys. Chem. B 2018, 122, 2867-2880) was applied to study antibody interactions. Two adjustable parameters related to the short-range interactions on the variable and constant regions were determined by fitting experimental data of 20 IgG1 monoclonal antibodies at 150 mg/mL. The root-mean-square deviation improved from 1 to 0.68, and the correlation coefficient improved from 0.63 to 0.87 compared to that of a previous model that assumed the short-range interactions were the same for all the beads. Our model is also able to calculate the viscosity over a wide range of concentrations without additional parameters. A tabulated viscosity based on our model is provided to facilitate antibody screening in early-stage design.

Differences in human IgG1 and IgG4 S228P monoclonal antibodies viscosity and self-interactions: Experimental assessment and computational predictions of domain interactions.

Human/humanized IgG4 antibodies have reduced effector function relative to IgG1 antibodies, which is desirable for certain therapeutic purposes. However, the developability and biophysical properties for IgG4 antibodies are not well understood. This work focuses on the head-to-head comparison of key biophysical properties, such as self-interaction and viscosity, for 14 human/humanized, and chimeric IgG1 and IgG4 S228P monoclonal antibody pairs that contain the identical variable regions. Experimental measurements showed that the IgG4 S228P antibodies have similar or higher self-interaction and viscosity than that of IgG1 antibodies in 20 mM sodium acetate, pH 5.5. We report sequence and structural drivers for the increased viscosity and self-interaction detected in IgG4 S228P antibodies through a combination of experimental data and computational models. Further, we applied and extended a previously established computational model for IgG1 antibodies to predict the self-interaction and viscosity behavior for each antibody pair, providing insight into the structural characteristics and differences of these two isotypes. Interestingly, we observed that the IgG4 S228P swapped variants, where the CH3 domain was swapped for that of an IgG1, showed reduced self-interaction behavior. These domain swapped IgG4 S228P molecules also showed reduced viscosity from experiment and coarse-grained simulations. We also observed that experimental diffusion interaction parameter (kD) values have a high correlation with computational diffusivity prediction for both IgG1 and IgG4 S228P isotypes.Abbreviations: AHc, constant region Hamaker constant; AHv, variable region Hamaker constant; CDRs, Complementarity-determining regions; CG, Coarse-grained model; CH1, Constant heavy chain 1; CH2 Constant heavy chain 2; CH3 Constant heavy chain 3; chgCH3 Effective charge on the CH3 region; CL Constant light chain; cP, Centipoise; DLS, Dynamic light scattering; Fab, Fragment antigen-binding; Fc, Fragment crystallizable; Fv, Variable domaing; (r) Radial distribution function; H1 CDR1 of Heavy Chain; H2 CDR2 of Heavy Chain; H3 CDR3 of Heavy Chain; HVI, High viscosity index; IgG1 human immunoglobulin of IgG1 subclass; IgG4 human immunoglobulin of IgG4 subclass; kD, Diffusion interaction parameter; L1 CDR1 of Light Chain; L2 CDR2 of Light Chain; L3 CDR3 of Light Chain; mAb, Monoclonal antibody; MD, Molecular dynamics; PPI Protein-protein interactions; SCM, Spatial charge map; UP-SEC, Ultra-high-performance size-exclusion chromatography; VH, Variable domain of Heavy Chain; VL, Variable domain of Light Chain.

Machine Learning Feature Selection for Predicting High Concentration Therapeutic Antibody Aggregation.

Protein aggregation can hinder the development, safety and efficacy of therapeutic antibody-based drugs. Developing a predictive model that evaluates aggregation behaviors during early stage development is therefore desirable. Machine learning is a widely used tool to train models that predict data with different attributes. However, most machine learning techniques require more data than is typically available in antibody development. In this work, we describe a rational feature selection framework to develop accurate models with a small number of features. We applied this framework to predict aggregation behaviors of 21 approved monospecific monoclonal antibodies at high concentration (150 mg/mL), yielding a correlation coefficient of 0.71 on validation tests with only two features using a linear model. The nearest neighbors and support vector regression models further improved the performance, which have correlation coefficients of 0.86 and 0.80, respectively. This framework can be extended to train other models that predict different physical properties.

Insights into Membrane Translocation of Protegrin Antimicrobial Peptides by Multistep Molecular Dynamics Simulations.

Protegrin-1 (PG-1) is a cationic arginine-rich antimicrobial peptide. It is widely accepted that PG-1 induces membrane disruption by forming pores that lead to cell death. However, the insertion mechanism for these highly cationic peptides into the hydrophobic membrane environment is still poorly understood at the molecular scale. It has previously been determined that the association of arginine guanidinium and lipid phosphate groups results in strong bidentate bonds that stabilize peptide-lipid complexes. It has also been suggested that arginine residues are able to drag phosphate groups as they insert inside the membrane to form a toroidal pore. However, whether bidentate bonds play a significant role in inducing a pore formation remains unclear. To investigate the role of bidentate complexes in PG-1 translocation, we conducted molecular dynamics simulations. Two computational electroporation methods were implemented to examine the translocation process. We found that PG-1 could insert into the membrane, provided the external electric potential is large enough to first induce a water column or a pore within the lipid bilayer membrane. We also found that the highly charged PG-1 is capable in itself of inducing molecular electroporation. Substitution of arginines with charge-equivalent lysines showed a markedly reduced tendency for insertion. This indicates that the guanidinium group likely facilitates PG-1 translocation. Potential of mean force calculations suggests that peptide insertion inside the hydrophobic environment of the membrane core is not favored. We found that formation of a water column or a pore might be a prerequisite for PG-1 translocation. We also found that PG-1 can stabilize the pore after insertion. We suggest that PG-1 could be a pore inducer and stabilizer. This work sheds some light on PG-1 translocation mechanisms at the molecular level. Methods presented in this study may be extended to other arginine-rich antimicrobial and cell-penetrating peptides.

Systematic Mutagenesis of Oncocin Reveals Enhanced Activity and Insights into the Mechanisms of Antimicrobial Activity.

Oncocin is a proline-rich antimicrobial peptide that inhibits protein synthesis by binding to the bacterial ribosome. In this work, the antimicrobial activity of oncocin was improved by systematic peptide mutagenesis and activity evaluation. We found that a pair of cationic substitutions (P4K and L7K/R) improves the activity by 2-4 fold (p<0.05) against multiple Gram-negative bacteria. An in vitro transcription / translation assay indicated that the increased activity was not because of stronger ribosome binding. Rather a cellular internalization assay revealed a higher internalization rate for the optimized analogs thereby suggesting a mechanism to increase potency. In addition, we found that the optimized peptides' benefit is dependent upon nutrient-depleted media conditions. The molecular design and characterization strategies have broad potential for development of antimicrobial peptides.

Free Energy Calculations of Microcin J25 Variants Binding to the FhuA Receptor.

Computer simulations were performed to study the antimicrobial peptide microcin J25 (MJ25), a 21-mer peptide with an unusual lasso structure and high activity against Gram-negative bacteria. MJ25 has intracellular targets. The initial step of MJ25 acquisition in bacterial cells is binding to the outer-membrane receptor FhuA. Molecular dynamics simulations were implemented to study the binding mechanism of MJ25 to FhuA and to search for important binding residues. The absolute binding free energy calculated from combined free energy perturbation and thermodynamic integration methods agrees well with experimental data. In addition, computational mutation analysis revealed that His5 is the key residue responsible for MJ25 and FhuA association. We found that the number of hydrogen bonds is essential for binding of MJ25 to FhuA. This atomistic, quantitative insight sheds light on the mechanism of action of MJ25 and may pave a path for designing active MJ25 analogues.

Coarse-grained modeling of protein second osmotic virial coefficients: sterics and short-ranged attractions.

A series of coarse-grained models, with different levels of structural resolution, were tested to calculate the steric contributions to protein osmotic second virial coefficients (B(22,S)) for proteins ranging from small single-domain molecules to large multidomain molecules, using the recently developed Mayer sampling method. B(22,S) was compared for different levels of coarse-graining: four-beads-per-amino-acid (4bAA), one-bead-per-amino-acid (1bAA), one-sphere-per-domain (1sD), and one-sphere-per-protein (1sP). Values for the 1bAA and 4bAA models were quantitatively indistinguishable for both spherical and nonspherical proteins, and the agreement with values from all-atom models improved with increasing protein size, making the CG approach attractive for large proteins of biotechnological interest. Interestingly, in the absence of detailed structural information, the hydrodynamic radius (R(h)) along with a simple 1sP approximation provided reasonably accurate values for B(22,S) for both globular and highly asymmetric protein structures, while other 1sP approximations gave poorer agreement; this helps to justify the currently empirical practice of estimating B(22,S) from R(h) for large proteins such as antibodies. The results also indicate that either 1bAA or 4bAA CG models may be good starting points for incorporating short-range attractions. Comparison of gD-crystallin B(22) values including both sterics and short-range attractions shows that 1bAA and 4bAA models give equivalent results when properly scaled to account for differences in the number of surface beads in the two CG descriptions. This provides a basis for future work that will also incorporate long-ranged electrostatic attractions and repulsions.