Computational models for studying physical instabilities in high concentration biotherapeutic formulations

ABSTRACT

Computational prediction of the behavior of concentrated protein solutions is particularly advantageous in early development stages of biotherapeutics when material availability is limited and a large set of formulation conditions needs to be explored. This review provides an overview of the different computational paradigms that have been successfully used in modeling undesirable physical behaviors of protein solutions with a particular emphasis on high-concentration drug formulations. This includes models ranging from all-atom simulations, coarse-grained representations to macro-scale mathematical descriptions used to study physical instability phenomena of protein solutions such as aggregation, elevated viscosity, and phase separation. These models are compared and summarized in the context of the physical processes and their underlying assumptions and limitations. A detailed analysis is also given for identifying protein interaction processes that are explicitly or implicitly considered in the different modeling approaches and particularly their relations to various formulation parameters. Lastly, many of the shortcomings of existing computational models are discussed, providing perspectives and possible directions toward an efficient computational framework for designing effective protein formulations.

KEYWORDS:

• Biotherapeutics
• drug formulation
• physical instabilities
• aggregation
• phase separation
• viscosity
• molecular modeling
• high concentration

Introduction

Protein-based therapeutics are at the frontier of development in pharmaceutical industry with a fast-growing global market, where monoclonal antibodies (mAbs) represent the largest class of biotherapeutics with more than 80 mAb drugs approved to date in the United States alone.^{1, 2} In addition to the challenges associated with their structural and functional design, protein solutions often exhibit physical instabilities such as aggregation and phase separation that arise from a complex interaction network among protein molecules with solution components. As current trends in biologics pipeline shift toward high concentration formulations, controlling protein instabilities is becoming more challenging. At elevated protein concentrations (>100 mg/mL), phenomena such as multi-body interactions and crowding exacerbate physical instabilities and might lead to other undesirable behaviors such as elevated viscosity and thermodynamic instabilities.³ While challenging, achieving stable high concentration protein formulations is necessary for both moving toward a patient-centric drug product and expanding the biologics drug market. As such, there is a need for rapidly advancing our understanding of the behavior of biotherapeutics at elevated protein concentrations.

Indeed, mitigating protein instabilities during the development of commercially viable biotherapeutics requires identifying optimal but phase-appropriate formulations. This entails exploring the space that governs the relations between formulation conditions and solution behavior. However, this formulation space is vast, where many parameters describing the solution conditions (e.g., protein concentration, pH, buffer, and excipients) are closely related to many protein properties such as hydrophobicity, charge distribution, morphology, and size. In fact, high concentration protein formulations constitute complex solutions, where formulation parameters are strongly interconnected to protein behavior such that a change in one parameter could cause contradictory effects on the relation between formulation and protein stability.^3,4 Moreover, due to limitations in material, time and resources availability during early-stage development (e.g., drug-candidate selection and preclinical development), a thorough experimental exploration of the formulation space becomes significantly challenging. In this regard, the implementation of fundamentally and statistically based computational models provides complementary tools for in-depth elucidation of the protein behavior, as well as for the subsequent identification of potentially relevant formulations. Specifically, these models can help design biologic drug formulations by: (1) constraining the formulation space to be experimentally investigated; (2) providing understanding of the underlying mechanisms for the different instability processes; and/or (3) identifying the mechanisms by which different solution components (or excipients) modulate protein behavior in the formulation.

Over the past two decades, an increasing number of studies focusing on the development and implementation of a variety of computational modeling tools have been reported. These studies have focused on understanding and/or predicting the behavior of protein solutions from either a biological or a biopharmaceutical standpoint.^5–12 As such, this review aims to provide a survey of the state-of-art of the in-silico application of a wide-range of computational models for effectively studying physical instabilities in protein solutions within the context of concentrated conditions. The models summarized here span various techniques and length-scales, ranging from atomistic simulations, coarse-grain representations, kinetic models, as well as novel approaches that combine resolutions from different molecular representations with other types of statistical and mathematical implementations. The review starts with an overview of the diverse classes of computational approaches that one commonly finds for evaluating the physical processes involved in destabilizing protein solutions. A particular emphasis is given on highlighting the range of length- and timescales that they can cover, as well as the underlying assumptions of each type of model. This overview aims to provide a summary of key physical considerations, practical and conceptual advantages, and missing components in the different classes of models. Thereafter, the different adverse thermodynamic and transport phenomena commonly affecting high-concentration protein formulations are revisited, exploring protein instabilities such as protein aggregation, phase separation and elevated solution viscosity. From a formulation perspective, protein–protein and protein–excipient interactions are the controlling knobs to modulate protein instabilities. Considering the landscape of mathematical models presented above, the manuscript continues with an overview of the proposed approaches used for evaluating protein interactions in both diluted and concentrated conditions. The review closes with a discussion about the perspectives and possible directions toward an efficient computational framework for designing effective protein formulations. This includes an analysis of the shortcomings of existing computational models in terms of computational cost, accessibility and inherent modeling limitations in capturing relevant experimental outcomes, as well as an examination of emerging multi-scale modeling approaches such as combining atomistic or coarse-grained models with machine learning or continuum models.

Types of protein models

Several computational models and tools have been developed in recent years, addressing many of the pharmaceutics-related protein stability problems such as protein self-association, protein aggregation, phase separation and elevated viscosity. These models and tools are very broad in terms of how proteins are represented and what “key” protein features are incorporated to study the different instability processes. Due to the multi-scale nature of protein stability, it is computationally prohibitive to use a single model to study both nanometer-scale problems (e.g., conformational changes, protein–protein interactions) and macroscopic issues (e.g., aggregation, phase-separation). As a result, the computational study of protein stability lends itself to a hierarchy of models (Figure 1), which can be classified as: 1) atomistic; 2) coarse-grain; and 3) continuum models. Although other types of modeling such as quantum mechanical representations and statistical approaches have been also used, this review primarily focuses on classical modeling studies applied to understanding and predicting protein stability phenomena, and problems particularly related to high-concentration protein formulations. Table 1 provides a summary of the different classes of models and their applicability to physical instability phenomena in protein solutions.

Table 1. Application of different model resolutions to protein instability processes

Resolution Level	Structure Prediction	Dynamics	Folding	Protein Interactions	Aggregation	Phase Separation	Rheology	Model Gaps
Atomistic	Yes	Yes	Yes^a	Yes^b	Yes^a	Yes^a	Yes^a	Large number of model parameters Lack benchmarking against pharmaceutically relevant conditions Limited to explicit simulations of small proteins and/or process with small characteristic timescales (<1 μs)
High-resolution coarse-grain	Yes	Yes	Yes	Yes	Yes^a	Yes^a	Yes^a	Large number of model parameters Lack benchmarking against pharmaceutically relevant conditions Limited resolution for side-chain specificity and dynamics Limited to timescales of ~1 ms
Low-resolution coarse-grain	No	Yes^c	Yes^c	Yes^d	Yes	Yes^a	Yes^a	Limited structural details and protein flexibility Non-existing relevant solvent/excipient models for this resolution level Applicable to length-scales on the order of ~100 nm Limited transferability between systems
Simplified coarse-grain	No	No	No	Yes^d	Yes	Yes	Yes	Only protein quaternary structure is preserved Non-existing solvent/excipient models for this resolution level No explicit information about structural features No transferability between proteins and formulations
Continuum	No	No	No	Yes^d	Yes	Yes	Yes	Minimal molecular-level details Limited to simplified protein representations (e.g., spheres, rods, ellipsoids) Applicable to macroscale processes Results limited to macroscopic properties

^aDue to computational cost, applications are generally limited to peptides and small proteins.

^bFor atomistic models, protein-protein/excipient interactions might be evaluated for any system; however, practical limitations make its use challenging for conditions leading to weak interactions.

^cAlthough some low-resolution CG models incorporate backbone flexibility, those degrees of freedom are generally treated from a mean-field approximation.

^dThe lack of structural resolution or explicit representation of one or more solution species limits the application of these models to capture specific interactions.

Figure 1. Representation of the hierarchy of computational protein models based on their level of resolution, using an IgG2 mAb (PDB: 1IGT) as an example. These types of models include: atomistic; high-resolution coarse-grain (based on the model from Bereau and Deserno¹³); low resolution coarse-grain from Blanco et al.;¹⁴ simplified coarse-grain using the 12-bead model from Calero-Rubio et al.¹⁵ and the 4- and 7-bead models from Blanco et al.;¹⁶ and a continuum model based on Wertheim’s theory adapted by Skar-Gislinge et al.¹⁷ The arrow indicates the direction in which the resolution-level increases for each model.

Atomistic models

Atomistic representations of both protein and all other species in the formulation (e.g., water, excipients) provide the most fundamental modeling approach to evaluate the behavior of protein solutions. In all-atom simulations, each atom is explicitly modeled as a single bead, and the protein solution behavior is characterized through the potential energy function. This function accounts for the bonded interactions (e.g., bond stretching, angle bending, torsions and improper angles), as well as intra- and inter-molecular non-bonded interactions (e.g., electrostatics, van der Waals, hydrogen bonding) for all molecular species. Consequently, the large number of different atomic forces yield an even larger number of parameters to capture the strength, range, and equilibrium energy for each type of interaction. The set of parameters or force fields for basic biomolecular systems (namely, proteins, nucleic acids, lipids, and biologically relevant ions) have been determined by several groups, where CHARMM¹⁸ and AMBER¹⁹ are among the most commonly used force fields. This class of models provide valuable insight into the early stages of different stability processes beyond experimental resolution capabilities. As such, they have been successfully applied to study various protein phenomena in crowded environments such as the kinetics and thermodynamics of protein conformational changes,²⁰ early stages of protein aggregation in peptides and small proteins,⁶ protein–protein/excipient interactions,²¹ among other processes.^22,23

However, there are some limitations of atomistic simulations in protein systems. Due to the significantly large number of particles or atoms that needs to be considered, in-silico experiments of protein solutions are generally constrained to processes occurring in timescales smaller than microseconds. While advances in enhanced sampling algorithms in tandem with the use of supercomputers have allowed for modeling large systems involving crowded environments and high protein concentrations,^24,25explicit simulations of phenomena of biopharmaceutical interest such as particle formation remains unreachable due to the long time- and length-scale for these stability processes.³ Similarly, due to typical parameterization strategies of atomistic models against experimental data, the inherent complexity of force fields leads to difficulties regarding transferability and relevance of these models between different types of biomolecular systems. Such issues become more evident when studying concentrated protein solutions, as the balance between protein–protein and protein-solvent/excipient interactions is not adequately captured.^26–28 Moreover, parametrization of force fields have been largely based on systems of biological interest rather than biopharmaceutical relevance, which poses a challenge when studying the effects of various formulation components such as polysorbates²⁹ and cryoprotectants.³⁰ Despite the substantial advances in recent years toward the improvement of force fields,^26–28 further efforts may be needed to fine-tune the different terms in the non-bonded interactions when simulating relevant protein formulations.

Coarse-grain models

Coarse-grain (CG) protein models have emerged as an alternative approach for studying protein stability problems, as they offer the potential for overcoming some of the inherent limitations of all-atom simulations.^12,13,31 In this class of models, the complexity of protein systems is simplified by grouping two or more atoms into a single particle (CG site), reducing the degrees of freedom in the system and expanding the range of time- and length-scales achievable by atomistic models.^7,9 Interactions between CG sites is however independent of the set of atoms used for mapping these sites, as they are developed to capture key physicochemical factors to successfully studying the phenomenon of interest. As such, the resulting mapping of a given set of atoms onto a CG site is neither unique nor arbitrary. This flexibility in coarse-graining proteins has yielded a wide range of CG models, from simplified models to highly detailed CG representations (Figure 1). Table 2 highlights selected models and their applications based on the hierarchy of coarse-graining resolutions.

Table 2. Selected examples of coarse-grain models and their applications

Model	Representation	Applications
High-Resolution Coarse-Grain
Bereau and Deserno^Citation13	3–4 beads per residue	Protein folding and peptide aggregation;^Citation13 Peptide-lipid interactions;^Citation32 Self-assembly of peptide copolymers.^Citation33,Citation34
MARTINI^Citation35	Up to 5 beads per residue	Dynamics of membrane proteins;^Citation36 Fibril growth in peptides;^Citation37 Self-assembly of peptides-DNA conjugates.^Citation38
OPEP^Citation39	Up to 6 beads per residue	Solution hydrodynamics;^Citation39,Citation40 Protein docking;^Citation41 Peptide fibrillation.^Citation42
PRIME^Citation43	4 beads per residue	Peptide aggregation in crowded environments;^Citation44 Aggregate polymorphism;^Citation45
PRIMO^Citation46	Up to 8 beads per residue	Membrane protein dynamics;^Citation46 Homology modeling and structure elucidation.^Citation47
Low-Resolution Coarse-Grain
Blanco et al.^Citation14	1 bead per residue	Protein self-assembly;^Citation48 Protein-protein interactions.^Citation49,Citation50
CABS^Citation51	Up to 4 beads per residue	Prediction of aggregation-prone regions as part of AGGRESCAN-3D;^Citation52 Protein-peptide docking;^Citation53 Protein folding.^Citation54
Kim and Hummer^Citation55	1 bead per residue	Structure refinement from SAXS data;^Citation56 Protein phase behavior;^Citation57 Protein-protein interactions.^Citation58
UNRES^Citation59	2 beads per residue	Structure prediction;^Citation60 Protein folding;^Citation61 Peptide aggregation.^Citation62
Simplified Coarse-Grain
Calero-Rubio et al.^Citation15	6 and 12 beads per mAb	Protein-protein interactions for mAbs at low- and high-concentration conditions.^Citation63–65
Chaudhri et al.^Citation66	12 and 26 bead per mAb	Solution dynamics and protein cluster formation in concentrated solutions.^Citation66–68
Wang et al.^Citation69	12 bead per mAb	Solution structure and Brownian dynamics of mAb solutions;^Citation69 Rheological properties.^Citation70
Vàcha and Frenkel^Citation71	1 capped-cylinder per peptide	Peptide self-assembly;^Citation71 Peptide aggregation at surfaces;^Citation72 Secondary nucleation during peptide fibrillation.^Citation73

In high-resolution CG models, residues are generally represented by four to seven CG sites. Most of these CG sites are assigned to the peptide backbone to preserve the dynamics of protein secondary structure, while the remaining sites represent the residue’s side chain to incorporate the identity of the protein sequence through amino acid-specific interactions. Examples of these models include the MARTINI model,³⁵ the OPEP model,³⁹ and the PRIMO model.⁴⁶ Because of the high level of structural details that these CG models provide, they have been successfully implemented for studying protein–protein interactions,⁷⁴ the mechanisms for protein folding and aggregation in small to medium size proteins and in crowded environments,^13,44 the self-assembly pathways of virus proteins,⁷⁵ as well as to facilitate the refinement of NMR and crystallographic structures.⁷⁶ High-resolution models represent a significant improvement over atomistic models in terms of computational cost, enabling the evaluation of protein process with characteristic timescales on the order of milliseconds.^9,77

By contrast, low-resolution CG models take the representation of proteins one step further by using only 1–3 CG sites per residue, where the specificity of different residue–residue interactions is maintained by explicitly incorporating the type and nature of each amino acid.^12,78 By reducing most of the degrees of freedom from the backbone, these models enable faster sampling of systems with multiple proteins or with characteristic length-scales on the order of ~100 nm.⁷⁹ In spite of the limited amount of structural details of the protein backbone, some of these CG models still account for the flexibility of the peptide bond (i.e., the distribution of distances, planar and torsional angles) from statistical analysis of the structural properties of known peptides and proteins.^59,77,78 For instance, the UNRES model,⁵⁹ a 2-CG site per residue representation, has been extensively used for studies of protein folding, structure prediction and the mechanism of protein fibrillation.^60–62 Other CG models in this category have been used to evaluate weak protein–protein interactions,^14,50 the effect of mutations on protein self-association,⁴⁸ the self-assembly mechanism of large protein complexes,^55,58 and the role of post-translational modifications on protein micro-phase separation.⁸⁰

On the other end of the spectrum of protein coarse-graining, there are simplified CG models, where fragments, or even entire proteins, are modeled by a single CG site. These models sacrifice sequence-level resolution to facilitate studying systems with as many as 10⁵ molecules.⁷¹ Nonetheless, elaborated force fields that include orientational-dependent interactions are developed for some of these type of models in order to capture key properties of the sequence heterogeneity. Moreover, their simplicity has enabled simulating systems of large proteins like mAbs without using state-of-the-art supercomputers. Indeed, a number of different simplified CG models for mAbs have been developed in recent years, where mAb representations span from 3 to 26 beads with different levels of intramolecular flexibility.^{15,16,66,81,82} These mAb models, as well as other simplified CG models for globular proteins, have been used for explicitly simulating macroscopic behavior of concentrated protein solutions such as crystallization, liquid–liquid phase separation, fibril formation, and their transport properties.^{70,73,83–85}

Similar to atomistic models, there are a number of different challenges in the development and implementation of coarse-grain models. These models are simplified representations of their all-atom counterparts; however, scaling between CG and atomistic models is not symmetric, as different CG sites within a given model might correspond to a different number of atoms with different chemical properties. Such asymmetry can result in both a non-uniform scaling of the system dynamics and a bias of the interactions between different CG sites, which may lead to difficulties in appropriately capturing the kinetics and thermodynamics of a given protein process.^7,86 Furthermore, CG models are generally “custom-made” toward studying a particular system of interest, where their parameterization is usually based on reproducing the system behavior at a given thermodynamic condition. As a consequence, it is often observed that these models cannot be transferred between different protein systems, or they are even unable to predict the behavior of the same system at different thermodynamic states.^82,87,88 Finally, most CG models are developed in an implicit-solvent framework, where the behavior of the solvent and any excipient in solution is averaged out and absorbed in the potential energy function for protein–protein interactions. If one were interested in studying the effects of excipients on the stability of protein solutions (e.g., during formulation screening studies), a different set of CG force-fields for each combination of excipients would be required to carry out such in-silico studies.

Continuum models

Continuum models have emerged as alternative tools for explicitly simulating protein processes with very large characteristic time- and length-scales (e.g., timescales larger than seconds and length-scales larger than micrometers), which are out of reach for atomistic and CG models.^89–91 These models aim to solve mechanistic thermodynamic, hydrodynamic and/or kinetic equations of a process of interest. However, unlike the other types of models discussed above, continuum models incorporate minimal molecular or structural details of the system to be studied. Instead, they generally rely on physicochemical properties of both the protein and solution such as diffusivities, surface tension, and various free-energies and kinetic rate constants, while molecular aspects such as protein–protein and protein-solvent/excipient interactions are often treated in a mean-field approximation.^92–94 This relative simplicity of the continuum models facilitates studying the behavior of protein solutions with significantly less computational resources than molecular models.

Among the different types of continuum models, those based in computational fluid dynamics (CFD) are arguably the most broadly used in the pharmaceutical industry, with various applications to different stages of upstream and downstream process development.^95–97 The success of CFD to study problems involving fluid flows has led to the development of a number of different models for studying the behavior of concentrated protein solutions under mechanical stresses such as shear forces in pre-filled syringes^92,98 and dense environments like subcutaneous tissue.⁹⁹ Similarly, several mechanistic kinetic models for protein aggregation have been developed to assess and predict the effect of different formulation parameters (e.g., pH, ionic strength, temperature) on the nucleation and growth of high molecular weight species. These aggregation kinetic models include mass-action kinetic approaches,¹⁰⁰ fixed-point approaches,¹⁰¹ and stochastic approaches.¹⁰² Other examples of continuum models include the application of theories based on statistical mechanics such as the Self-Consistent Ornstein-Zernike Approximation, Kirkwood-Buff solution theory, and Wertheim’s perturbation theory for evaluating high-concentration phenomena such as protein self-assembly,¹⁰³ liquid–liquid phase separation,¹⁰⁴ protein interactions,^105,106 and solution rheology.^107,108

Overall, continuum models are advantageous in connecting in-silico analysis of the behavior of protein solution to experimental development of biologic drugs. These models focus on evaluating stability phenomena on a scale comparable to that of most experimental techniques. In fact, they are often used for fitting experimental data to augment the information obtained from different assays, as well as to predict the behavior of protein formulations at different solution conditions.^85,107,109 However, special considerations are needed when implementing these types of modeling. While the underlying theories that constitutes most continuum models are rigorous, multiple assumptions are required to adapt these theories for modeling complex systems such as protein solutions. Due to these assumptions, the accuracy and relevance of the models may only hold for a small subset of solution conditions, which might result in misinterpreting the behavior of the solution beyond the conditions directly compared against experiments. For instance, a common simplification in statistical mechanics theories is to neglect high-order dependence of protein concentration.¹⁰⁴ Although such assumption is valid at diluted conditions, it fails to correctly describe the behavior of concentrated solutions due to factors like multi-body interactions and long-range spatial correlations.¹¹⁰ Likewise, because of the large simplification in the molecular resolution of the system, the individual effects from different physical factors (e.g., different types of interactions, molecular anisotropy and heterogeneity) are reduced to a few parameters in the continuum models, leading to potential challenges in analyzing the results from the models.¹⁷ The balance between these physical factors is a function of the local protein environment, and thus one should expect that the reduced model parameters are not constant but change with the solution conditions.⁸² Additionally, the roles that the different physical factors play on the behavior of the solution are not independent from each other, which might lead to multiple combinations of the model parameters to describe the behavior of the protein system at a given solution condition.¹⁰⁰ Because of these challenges, the degree of predictability and transferability of the continuum models to different protein systems and formulations is very limited.

High-concentration physical instabilities

Protein self-association and aggregation

Among the different physical stability issues affecting proteins, aggregation is indisputably the most prevalent problem during the development of biotherapeutics. The formation of high molecular weight species and particulate matter can reduce the efficacy and affect the appearance of the drug product, in addition to potentially lead to unwanted immunogenicity.³ Predicting and minimizing protein aggregation is therefore an important drug development goal. However, this remains a challenging task, as aggregation is an ubiquitous process to any protein, where the rate for aggregate formation is very sensitive to both protein structure and solution conditions.^111,112 Moreover, aggregation is a complex multistep process governed by intra- and intermolecular interactions, with characteristic time and length scales that span many orders of magnitude.^100,111

Mechanistically, protein aggregation occurs through a series of both reversible and irreversible stages^100,111 (Figure 2), which include: conformational change of the protein monomer to form an aggregation-prone or reactive species (stage I); nucleation via protein self-association (stage II); the formation of the smallest irreversible aggregate species (stage IIIa) or a homogeneous phase separation when native self-association occurs (stage IIIb); aggregate growth via monomer addition (stage IV), aggregate–aggregate coalescence (stage V), and aggregate fragmentation (stage VI); and the phase separation or precipitation of the high molecular weight species (stage VII). Depending on the specific protein system, the relevance and rates of these stages can differ, leading to a variety of aggregation mechanisms. In this regard, in-silico models provide invaluable tools to gain insight into the critical factors affecting the different aggregation stages. Although most of the existing computational models cannot capture the complete range of relevant time and length scales, different types of modeling approaches have been developed and implemented to independently assess these stages, making predictions of the aggregation propensity and long-term stability of a given protein formulation.^86,113,114

Figure 2. Schematic representation of the generalized protein aggregation mechanism for multidomain proteins such as mAbs. The stages shown in the diagram correspond to either effectively reversible steps (double arrows) or irreversible steps (single arrows). Protein oligomerization can occur through self-association of the native monomer (N) or a (partially) unfolded reactive species R. The mechanism also considers the case that N self-associates to a critical size (N_X) to nucleate a homogeneous phase separation (e.g., liquid-liquid separation or crystallization).

Based on the generalized aggregation reaction in Figure 2, different (continuum) aggregation kinetic models have been established.^{93,100,102,109,115,116} All these models express the extent of aggregation via mass balance equations for the elementary reactions, but they differ by how the aggregate size distribution is treated. One of the most comprehensive aggregation models is that developed by the Roberts’ group, which follows a Lumry-Eyring Nucleated-Polymerization approach (LENP).¹⁰⁰ This model incorporates most of the aggregation stages outlined above, while it explicitly captures a broad range of aggregate species to yield a discrete aggregate size distribution. The LENP model has been successfully applied to characterize the effect of solution conditions on the aggregation mechanism of globular proteins¹¹⁷ and antibodies.¹¹⁸ However, implementation of the LENP model is somewhat limited to diluted conditions, as relevant factors to the mass transport in concentrated protein solutions (e.g., crowding, viscosity effects, ion binding) are not incorporated in the model. Alternative approaches combine the aggregation mass balance equations with the Smoluchowski coagulation equations, which allows for coupling the elementary rate constants to the diffusion and fractal dimension of the aggregated species.^{114,116,119,120} This type of kinetic model has been shown to fit experimental data for mAbs reasonably well over a wide range of solution conditions and protein concentrations as high as 60 mg/mL.¹¹⁶ Nonetheless, in a modified Smoluchowski-based model that incorporates concentration-dependent viscosity, Nicoud et al.¹⁰⁹ found for a concentrated mAb solution that large discrepancies between the aggregation model and experiments are obtained when strong anisotropic protein interactions and/or native aggregation are suspected. On the other hand, other aggregation kinetic models have focused on solving the mass balance equations in terms of probability density functions¹²¹ and generating functions^102,115,122 to recover the aggregate size distribution. These models have successfully captured the time-evolution at various protein concentrations^101,102 and at crowded conditions.¹²¹ Moreover, some of these models, as those developed by Knowles and collaborators,^101,122,123 yield a simplified analytical solution to the set of reaction equations, which reduces the computational burden and facilitates their implementation for fitting experimental data. However, most, if not all, of these latter kinetic models have been exclusively implemented to study fibrillation in peptides and small proteins. As such, they generally simplify aggregate growth to occur only via monomer addition, while other stages commonly observed during amyloid formation (e.g., secondary nucleation, auto-catalysis, and stochastic kinetics) are incorporated.

An alternative computational approach is also commonly used to evaluate and predict aggregation propensity and kinetics, which is based on information derived from single-molecule features. Different experimental mutagenesis studies on amyloidogenic proteins have shown that aggregation is correlated to physicochemical features in protein structure such as hydrophobicity, charge, β-strand propensity, and surface arrangement of amino acids.^124,125 These observations initially spurred the development of a first-generation of aggregation predictors based solely on protein sequence, which identify aggregation-prone regions (APRs) from pattern matching and heuristic equations validated against experimental databases of hundreds of amyloid-forming short peptides (Table 3). Among these predictors, there are methods such AGGRESCAN,¹²⁶ TANGO,¹³⁰ PASTA,¹²⁹ and WALTZ.¹³¹ The interested reader is referred to recent reviews on this specific topic.^86,139 These sequence-based algorithms provide fast and computationally inexpensive tools to identify APRs and rank proteins based in their intrinsic aggregation propensity, and they have been used for guiding modification of mAb sequences to improve their stability.¹⁴⁰ However, these predictors present several limitations for a broader applicability in the biopharmaceutical industry. Firstly, most of these sequence-based tools predict APRs based on the aggregation mechanism of peptides and amyloidogenic proteins, which may not be representative of the aggregation pathway of larger proteins. Secondly, available experimental data for building and validating these algorithms is limited to a few solution conditions (mainly, physiological conditions), and therefore more data is required to expand their usefulness for predicting APRs in biotherapeutics over relevant formulations (e.g., larger pH range, different buffers and excipients). Lastly, the inherent assumptions of these sequence-based algorithms prevent them from capturing the role that the three-dimensional structure plays on aggregation. In fact, when applied to large proteins such as mAbs, these predictors often lead to a larger number of false positive and negative results, identifying regions that are not solvent exposed and failing to capture APRs involving residues from non-continuous sequence fragments.^8,141

Table 3. Selected examples of aggregation-propensity predictors

Model	Underlying Approach	References
Sequence-Based APR Predictors
AGGRESCAN	Aggregation propensity scale per residue based on in-vivo experiments of amyloidogenic proteins	Conchillo-Solé et al.^Citation126
AmyloGram	Aggregation propensity derived from 21 amino acid physicochemical properties such as size, hydrophobicity, polarity, secondary structure propensity and contact propensity	Burdukiewicz et al.^Citation127
PAGE	Aggregation propensity derived from physicochemical properties including aromaticity, β propensity, polarity, charge and solubility	Tartaglia et al.^Citation128
PASTA 2.0	Aggregation propensity derived from hydrogen-bonding energy functions for non-bonded residue pairs in beta-strand structures	Walsh et al.^Citation129
TANGO	β-sheet formation propensity based on empirical and statistical correlations derived from energy functions	Fernandez-Escamilla et al.^Citation130
WALTZ	Aggregation scale per residue derived from the aggregation propensity of a set of hexapeptides	Maurer-Stroh et al.^Citation131
Zyggregator	Aggregation prediction based on heuristic correlation accounting for hydrophobicity, secondary structure propensity, net charge, and presence of Gatekeeper residues	Tartaglia and Vendruscolo^Citation132
Structure-Based APR Predictors
AGGRESCAN-3D	Aggregation propensity derived from combining AGGRESCAN score with calculated exposed surface area from simulations with CABS CG model^Citation53	Kuriata et al.^Citation52
AggScore	Solvophobic patches calculated from atomic partial charges and logP values	Sankar et al.^Citation133
CamSol	Aggregation prediction derived from combining the score from Zyggregator with solubility calculations based on protein structure	Sormanni et al.^Citation134
Developability Index	SAP values together with residue charge calculations derived from predicted residue pKa values	Lauer et al.^Citation135
SAP	Residue hydrophobicity, solvent accessible area derived from all-atom Molecular Dynamics simulations	Chennamsetty et al.^Citation136
SolubiS	Aggregation propensity based on APR prediction from TANGO and conformational stability calculations from the FOLDX force field^Citation137	Van Durme et al.^Citation138

To overcome some of the limitations of the sequence-based tools, a second generation of algorithms for predicting APRs have emerged, which explicitly accounts for the folded structure of proteins. That is, rather than relying solely on the protein sequence, the protein structure (e.g., from crystallographic or homology models) is leveraged to assess the likelihood that a given APR might be involved at the interface for aggregate formation. Some of these structure-based predictors, such as SolubiS,¹³⁸ CamSol,¹³⁴ and AGGRESCAN-3D,⁵² combine the predictions from sequence-based tools with calculations of conformational stability and residue energy to weight APRs based on their tendency to be solvent-exposed or to interact with other fragments. Other predictors like SAP,¹³⁶ Developability Index,¹³⁵ and AggScore,¹³³ predominantly use the 3D protein structure to assess the residue solvent accessibility and/or partial charges (e.g., via short Molecular Dynamics simulations), which are then correlated to aggregation propensity. As a result, these structure-based predictors do not only outperform their sequence-based counterparts, but they have also been found useful for screening and re-engineering biotherapeutics candidates to reduce their aggregation propensity.⁸ Moreover, Wolf et al.¹⁴² found that the results obtained from some of these in-silico algorithms for a series of mAbs are well correlated with several experimental techniques used to evaluate their early-stage developability. In other studies, Trout and collaborators have combined the calculation from SAP with Molecular Dynamics simulations of protein-excipient interactions to gain insight into how aggregation propensity and viscosity of mAb formulations are affected by excipients such as carbohydrates³⁰ and ionic species.¹⁴³ Despite the significant improvements of the structure-based predictors, there are a few downsides that need to be considered for their implementation. Calculations of the different structural properties are carried out from fluctuations around the native protein structure, and thus they inherently bias the results toward APRs involved in a native-aggregation pathway. Likewise, these algorithms are based on single-molecule simulations, where the general assumption is that each APR is independent of each other in terms of their contribution to protein aggregation. As a consequence, they neglect synergistic effects that may arise from the proximity of two or more APRs.⁸ The single-molecule calculations also make impossible to distinguish protein concentration effects on aggregation such as the role of multibody interactions and excluded volume.¹¹¹ Nevertheless, as more experimental data on aggregation of biotherapeutics becomes available, these limitations might be overcome by improvements in the underlying heuristics to correlate structural properties to aggregation propensity. Promising efforts in this direction have been recently seen by Lai et al.,¹⁴⁴ where different structural properties, including those used in SAP, were fed into a machine learning algorithm to predict the aggregation rate of 21 mAbs at high-concentration.

The use of molecular models has not been limited to investigating the molecular properties related to aggregation. Both atomistic and CG models have been widely used to understand the dynamics and thermodynamics of various aspects of the aggregation mechanism.^6,78,113 The development of enhanced sampling methods such as Replica-Exchange Molecular Dynamics, Metadynamics, Umbrella Sampling, and Markov modeling, among others, have facilitated expanding the use of these different types of protein models to study critical steps from the generalized aggregation mechanism (Figure 2). Reviews that provide a comprehensive summary of advanced computational methods used for studying protein aggregation and other processes are available.^10,24 In the case of all-atom simulations, studies on protein aggregation are often focused on small protein fragments from known APRs. In-silico studies of homopeptides (e.g., poly-alanine, poly-valine, poly-glycine) and several fragments from amyloidogenic proteins have enabled our understanding, at a molecular level, of the early-stage mechanism of fibril formation.^145–147 These studies have shown that peptides initially collapse into a partially ordered oligomeric state up to a critical size nucleus of 6–8 strands, to then evolve into ordered β-sheet structures. Luiken and Bolhuis¹⁴⁸ showed that this nucleation process can change from a one-step to a two-step nucleation mechanism as the peptide hydrophobicity increases. Atomistic models have also provided insight into the growth stages of fibrillation.^6,31,113 These models have highlighted the role that conformational fluctuations play during elongation^149,150 and secondary nucleation,¹⁵¹ as well as enabled the characterization of key residue–residue and residue–water interactions that govern the stability and fragmentation of fibrils from different amyloid-prone peptides.^152,153 Note that the implementation of atomistic models have mainly focused on amyloid formation, whereas their application to the phenomenon of aggregation in biotherapeutics have been directed to characterize the local dynamics of antibody fragments to help explain experimental observations of different aggregation behaviors.¹⁵⁴

On the other hand, computational studies of the aggregation mechanism of larger proteins and/or longer length-scales have been possible via coarse-grain models from all different resolution-levels.^7,31,155 Different simplified CG representations have been developed to elucidate how the interplay between soluble and aggregation-reactive conformations affect the nucleation and elongation stages of fibrillation.^71,155,156 In an interesting report by Ŝarić et al.⁷³ using one of these simplified CG models, the authors found that monomers can spontaneously aggregate without a nucleation step at high protein concentrations, while the formation of a small oligomeric nucleus is a prerequisite for aggregation at low protein concentrations. Likewise, low-resolution CG models have led to new insights into protein self-association and aggregate growth.⁷⁸ Phenomenological models such as those from Shea’s and Caflisch’s groups have been used to describe the effect of hydrophobic and charge residues on the nucleation stage,¹⁵⁷ as well as the dependence of the different pathways for fibril growth on protein conformation.^158,159 Other low-resolution CG models based on rigid representations of proteins have been used to understand the role of surface residues and solution conditions of protein self-association.^49,50,62 For instance, Blanco et al.¹⁴ used a 1-bead per residue protein model for γD-crystallin to evaluate the effect of ionic strength on modulating preferential protein orientations during self-association, enabling the identification of key mutations to reduce the aggregation rate.⁴⁸ Notably, protein aggregation studies using different high-resolution CG models have arrived at similar conclusions regarding the mechanisms of nucleation and aggregate growth (albeit using smaller fragments of the same proteins).^40,160,161

Liquid-liquid phase equilibrium

Another concerning physical stability issue in high-concentration protein formulations is their potential to become opalescent and undergo a liquid–liquid phase separation (LLPS) process during refrigerated conditions, where the solution separates into protein-rich and protein-poor phases.^11,162,163 This phenomenon can both affect the esthetics of the drug product and trigger other stability issues.¹⁶³ While LLPS is typically a reversible process (e.g., increasing temperature brings the solution back to a homogeneous, single phase³) the partitioning of the different solution components (particularly, ionic species) between both phases might trigger other phenomena. The imbalance of the buffer/excipient species between the equilibrium phases might shift the pH and ionic strength toward unfavorable conditions and result in protein unfolding, irreversible aggregation or protein precipitation.¹⁶³ This problem is not exclusive to the biopharmaceutical industry, as LLPS also occurs in living cells and is related to mechanisms of intracellular organization and various diseases.^164,165 As such, understanding the phase separation of protein solutions remains an active research area in many disciplines from both experimental and computational standpoints.

Experimental studies on several globular proteins such as lysozyme¹⁶⁶ and γ-crystallins¹⁶⁷ have provided a comprehensive picture of the phase behavior of these proteins, with prominent common features such as a metastable LLPS with respect to crystallization and the formation of an arrested or gel-like state above the critical concentration (Figure 3). The resemblance of the phase diagram of globular proteins with that of short-ranged attractive colloidal particles spurred most of the earliest computational studies in LLPS for proteins.¹⁶⁸ In fact, simplified CG models consisting of hard-spheres with an isotropic short-range attraction are able to qualitatively reproduce both the metastability of the liquid–liquid transition and the high-concentration arrested state provided that the range of the intermolecular interactions is sufficiently small (e.g., ~1/8 of the protein diameter).¹⁶⁹ Moreover, these simplified colloidal protein models have led to the discovery of an “extended law of corresponding states” for LLPS in globular proteins.^169,170 That is, the LLPS of globular proteins collapses into a master curve when representing the binodal curve in terms of the strength of the attractive interactions (via the osmotic second viral coefficient, $B_{22}$) and an effective protein volume that accounts for the screened protein charge. However, these isotropic models are unable to quantitatively capture the shape of the phase boundaries and the concentration for crystallization, and instead CG models with patchy or directional interactions are required for better quantitative agreement with experimental data.^11,162 Numerous groups have not only demonstrated that these patchy models yield broader binodal curves like those of globular proteins,^171–173 but they also provide insight into the variety of space groups on protein crystals.⁸⁴ Moreover, it has been shown through these patchy models that LLPS can be either suppressed or triggered by protein oligomerization based on the extent of aggregation.^83,171 The computational study of LLPS has not been limited to simplified CG models, as other CG models with different levels of resolution, as well as continuum models, have also been developed for evaluating LLPS in globular and disordered proteins.^104,165,174 In a very interesting approach, Wertheim’s thermodynamic perturbation theory was adapted to capture directional interactions of globular proteins,¹⁰⁴ and it was even extended to study protein-phase separation in the presence of different buffers¹⁷⁵ and salts.¹⁷⁶

Figure 3. Generic phase diagram for globular proteins adapted from Muschol and Rosenberger.¹⁶⁶ The regions below the solubility curve (i.e., the gel and liquid-liquid coexistence regions) are metastable with respect to crystallization. The liquid-liquid coexistence region, bounded by the binodal curve, corresponds to the thermodynamic state where the solution separates into protein-rich and protein-poor phases. The gelation curve indicates the boundary for the formation of an arrested state. For any protein, the relative position between the solubility, binodal and gelation curves depends on both the protein sequence and solution conditions. Redrawn from Ref. 166.

In-silico studies of LLPS have also been applied to proteins of pharmaceutical interest such as monoclonal antibodies.^177,178 The liquid–liquid binodal curve of mAbs significantly differs from that of globular proteins, as the critical point of antibody solutions typically occurs at lower temperatures and concentrations, while the binodal curve is broader.^11,163 Sun et al.¹⁷⁷ used simplified CG models representing mAbs as flexible molecules composed of 3 to 7 CG-sites to demonstrate that both the ‘Y’ shape and the flexibility of the hinge region contribute to the asymmetrical shape of the binodal, but they have minimal effect on determining both the critical concentration and elevated density of the protein-rich phase. Instead, it was found that the inner subdomains (i.e., C_H1, C_H2 and C_L) need to be net repulsive or less attractive in comparison to the other subdomains in order to obtain a phase coexistence curve quantitatively comparable to experiments. Hinge flexibility, on the other hand, still plays an important role in facilitating quaternary structural rearrangement to achieve a more compact but stable solution structure at very high mAb concentrations.¹⁶ Recently, Vlachy and collaborators¹⁷⁸ arrived at similar conclusions by extending Wertheim’s theory for a 7-CG-site mAb model. These authors found that the critical temperature and concentration are sensitive to the imbalance of the interactions involving the C_H3 and variable fragment (F_V), but they are marginally affected by the actual strength of the intermolecular interactions. This latter continuum model has been implemented to semi-quantitatively reproduce experimental LLPS data for two different mAbs,¹⁷⁸ as well as to evaluate the effect of polymers¹⁷⁹ and bulky agents¹⁸⁰ on the phase separation of antibodies.

The development and implementation of both simplified CG and continuum models has driven much of our understanding of the phase behavior of proteins. These models have not only provided insight about the relevance of anisotropic or patchy interactions on the LLPS of proteins, but they have also allowed us to identify simple guidelines to qualitatively identify conditions that lead to phase separation (e.g., based on $B_{22}$ for globular proteins or the imbalance of interactions between mAb fragments). Nonetheless, further research is still needed to streamline the use of modeling for robustly screening protein drug candidates or drug product candidates against phase behavior. Firstly, while patchy models represent the state-of-the-art for quantifying LLPS, the definition of a protein “patch” remains loose. There is no comprehensive rationale to select the degree of anisotropy to connect the CG patchy model to the protein structure or sequence. Commonly, patches are placed in either a random or symmetric fashion, and experimental data is fit based on the number, size and interacting strength of these patches.¹⁰⁶ Although this approach can be effective, it does not give insight about specific molecular features related to LLPS for a given protein or class of proteins. Different efforts have been made to relate surface features or sequence fragments to interacting patches,^181,182 but it is unclear what the relevant characteristics are for these patches to influence LLPS and other instability phenomena. In this regard, the use of higher resolution CG models might provide an alternative approach to overcome these issues, as they can provide information about residue-level interactions related to LLPS. In fact, such models have been recently used to assess the LLPS of intrinsically disordered proteins.^164,165 Another outstanding challenge comes from the inherent modeling limitations for evaluating phase separation processes. From a computational perspective, identifying and characterizing phase coexistence curves are some of the most expensive and intensive modeling tasks, as they require sampling over millions of configurations and/or very large system volumes to overcome the generally large thermodynamic barriers between the phases at equilibrium. Moreover, these simulations typically need to be carried out over several thermodynamic states (e.g., multiple sets of temperature and pressure) in order to reconstruct the binodal curves. As such, novel and clever approaches need to be developed or adapted to reduce the computational burden for assessing LLPS. Some recent methodologies based on multi-scaling,¹⁸³ Widom insertion,¹⁸⁴ and thermodynamic extrapolation¹⁸⁵ might provide a path forward in this regard. Last but not least, there is a lack of experimental data for phase behavior of proteins, and in particular for pharmaceutically relevant biologics. LLPS data have only been reported for a few proteins, which represents a challenge for identify generalized guidelines to predict potential problems for a given drug candidate with respect to phase behavior.

Transport properties of high-concentration solution

Solution viscosity is a critical attribute for the development of high-concentration protein formulations, as an elevated viscosity (>30 cP) can significantly impact the pressure and flow in various unit operations such as filtration, ultrafiltration-diafiltration, and filling,^95,96 as well as limit the development of devices for drug administration (e.g., auto-injectors and pre-filled syringes).¹⁸⁶ From a physicochemical standpoint, it has been proposed that the presence of transient protein clusters is the root cause of high viscosity in protein solutions, which is in turn driven by protein–protein interactions.^187,188 As such, understanding the relationships between interactions, solution structure and solution rheology is key for establishing appropriate formulation strategies to achieve a suitable viscosity in high-concentration biotherapeutics. In this regard, computational protein models have greatly contributed to our current knowledge of the molecular origins of viscosity behavior in protein solutions.

Early computational work focused on investigating the viscosity problem using globular proteins as model systems, where well-established (continuum) colloidal models for spherical particles such as the mode-coupling theory¹⁸⁹ (MCT) could provide insight into the relation between protein interactions and solution viscosity. Several groups have studied the rheological behavior of solutions of bovine serum albumin (BSA) up to 300 mg/mL at different buffer conditions and salt concentrations.^190–192 Interestingly, colloidal models based on particles with hard-sphere repulsion and long-range electrostatic interactions were able to fit reasonably well the concentration-dependent viscosity of BSA, indicating that long-range repulsive interactions govern the behavior of viscosity.¹⁹¹ A similar conclusion was also found by Foffi et al.¹⁹³ for α-crystallin, where screening of electrostatic repulsions allows the solution viscosity to be captured by a simple polydisperse hard-sphere model and by MCT. Despite the success of these repulsive colloidal models, many proteins interact through a combination of short-range attractions and long-range repulsions. Different studies with lysozyme have demonstrated that the competition between these types of interactions facilitates the formation of protein clusters through an intermediate-range order, where screening electrostatics yields a higher solution viscosity.^107,110,194 In such cases, previously used colloidal models such as MCT fail to capture the viscosity of concentrated lysozyme solutions.¹¹⁰ In an interesting recent report, von Bülow et al.¹⁹⁵ used all-atom Molecular Dynamics simulations to evaluate the concentration dependence of viscosity for four small globular proteins (≤ 14 kDa) for protein concentrations as high as 200 mg/mL. The authors found that protein crowding strongly affects both translational and rotational diffusions, but the slowdown on the rotational diffusivity is mainly related to the formation of weak, dynamic protein clusters with dissociation constants of ~20 mM. Indeed, they derived a heuristic model that well reproduces experimental diffusivities based on the mean cluster size and viscosity calculated from the simulations.

When modeling the viscosity of pharmaceutically relevant proteins such as mAbs, the situation is more dire than in the case of globular proteins. The multi-domain nature and anisotropic shape of mAbs make the implementation of existing colloidal spherical models difficult for predicting solution viscosity. Moreover, under identical changes in formulation conditions, some antibodies show opposite trends in viscosity behavior.^4,188 These observations have suggested that the elevated viscosity of mAb formulations is driven by local sequence and structural features rather than net colloidal effects. As such, a large portion of the computational research carried out in this area has focused on the implementation of all-atom models to correlate local molecular descriptors (e.g., protein charge, solvent accessible area, dipole moment) with experimental viscosity data.⁸ In one of the earliest works with mAbs, Li et al.¹⁹⁶ evaluated a series of molecular descriptors for 11 homology models of mAbs based on electrostatic and solvophobic properties. The authors correlated these descriptors with viscosity measurements of the antibodies under equivalent solution conditions, finding that the viscosity behavior for a given mAb isotype is correlated to descriptors associated with the F_V domain, such as charge, pI, zeta potential, and aggregation propensity. Likewise, Tomar et al.¹⁹⁷ used experimental measurements of 16 different mAbs to develop a computational scheme to predict their concentration-dependent viscosity based on the electrostatic and solvophobic properties of both the F_V and the full mAb structure. In addition to the electrostatic properties of the F_V, this latter work identified that both the hydrophobic surface area of the mAb and the charge of hinge region are important to the predictability of viscosity. In agreement with these observations, Sharma et al.¹⁴¹ used 14 mAbs to correlate experimental viscosity values with in-silico molecular descriptors, including properties calculated from Molecular Dynamics simulations. The authors concluded that the solution viscosity was correlated with the hydrophobicity and charge dipole of the F_V region. More recently, Lai et al.¹⁴⁴ used 27 mAbs approved by the Food and Drug Administration to expand this type of modeling analysis, combining molecular descriptors obtained from atomistic models with machine learning feature selection. Both the net charge of the full mAb and the number of solvophobic residues in the F_V region were highlighted by the machine learning algorithm as key features for the viscosity behavior.

An alternative computational approach has focused on using simplified CG models for mAbs to study the relation between protein interactions and solution structure, seeking to correlate changes in the spatial correlations of mAbs with experimental measurements of viscosity. Chaudhri et al.^66,67 used two different simplified CG resolution levels to study the solution structure for two antibodies, which differ by only a few mutations in the complementarity-determining regions (CDRs) but exhibit different concentration-dependent viscosity. One of the CG models represents mAbs with one bead per subdomain (i.e., a 12-CG-site mAb representation), while the other model adds further resolution at the CDR and hinge regions by using 26 CG-sites to represent mAbs. Based on calculations of the potential of mean force at different protein concentrations via Molecular Dynamics simulations, it was found that complementary electrostatic interactions involving both the antigen-binding fragment (Fab) and F_C leads to the formation of antibody clusters, networks, and higher-order structures, which correlates with the viscosity behavior of both mAbs. Interestingly, the use of a higher level of coarse-graining did not provide much further insight into the solution structure of these mAb systems, as the 12-CG-site model was sufficient to determine the underlying cause of the concentration-dependent viscosity. A similar conclusion was also obtained when further implementing the same modeling approach on 4 additional mAbs.⁶⁸ Wang et al.⁶⁹ combined the 12-CG-site representation with Brownian dynamics simulations to semi-quantitatively reproduce the experimental transport properties (e.g., self-diffusivity, structure factor and viscosity) of the same mAbs used by Chaudhri et al.⁶⁶ The authors identified as the cause of changes in transport properties the formation of weakly interacting protein clusters rather than dense or strong networks. Nonetheless, it was required to impose nonphysical constraints to maintain rigidity in the protein cluster in order to reproduce the high-concentration viscosity behavior of one of the mAbs. This model was recently improved by Lai et al.⁷⁰ by incorporating anisotropy into the short-range interactions between the constant and variable regions of the mAbs instead of an uniform van der Waals interaction term for all the CG-sites. The improved 12-CG-site model was evaluated against experimental viscosity measurements of 27 antibodies, where model parameterization was directly coupled to a previous machine learning approach.¹⁴⁴

In a different set of works, Dear et al.⁸⁵ and Chowdhury et al.⁸⁸ used a similar 12-CG-site mAb model in combination with small-angle X-ray scattering (SAXS) experiments to predict the concentration-dependent viscosity of two mAbs over a broad range of formulations. SAXS data was used to determine whether anisotropic protein interactions were relevant at a given solution condition, as well as to parameterize the CG force field. The resulting CG model was used for calculating the cluster size distribution of the solution, which in turn was used to reproduce the solution viscosity via an empirical equation.⁸⁸ This approach predicted reasonably well the changes in viscosity with respect to protein concentrations for both mAbs in most of the tested formulations, as well as for a polyclonal IgG;¹⁹⁸ however, it failed to capture these changes when protein clustering is driven by strong anisotropic interactions. More recently, Izadi et al.⁸¹ developed a 10-CG-site mAb model, where the F_C domain is represented by only 2 CG-sites. Unlike all previous CG models, parameterization of the CG force field was based solely on data from atomistic simulations to incorporate the multipole moments of the charge distribution. The resulting model was compared against experimental data from 16 antibodies in two different formulations (low and high ionic strength), and it was qualitatively correlated with transport properties such as the diffusion interaction parameter and viscosity. Other recent models based on Wertheim’s theory have arrived at similar conclusions, showing promising results for semi-quantitatively capturing the relation between anisotropic interactions (between F_V and C_H3 domains) and solution viscosity.^17,199

Although both atomistic and CG models have been fundamental to expand our understanding of the role that electrostatic and anisotropic interactions play in the concentration-dependent viscosity of mAbs and other proteins, these models have some limitations, and as such, present opportunities for improvement. Benchmarking and validation of the different models have been limited to no more than 27 different biotherapeutics and only a few different solution conditions. A larger and more diversified set of molecules and formulations is needed to extend the validity, accuracy, and robustness of any predictive scheme. Additionally, the connection between protein interactions and viscosity behavior across the majority of the models have relied on either heuristic equations or statistical analysis. While these approaches allow us to establish correlations between the molecular and macroscopic behavior of the protein solution, they do not identify the underlying cause. Development of further phenomenological theories and models is still needed to elucidate the molecular origins for the concentration-dependent viscosity of a given biologic formulation. Moreover, most of existing models for assessing the viscosity behavior of protein solutions uses implicit solvent approximations. This approach facilitates reducing computational cost, but it significantly limits our ability to rationally identify the type of viscosity-modifier excipients that can be used for a given drug product formulation. While the work by von Bülow et al.¹⁹⁵ has demonstrated that it is possible to use all-atom modeling for explicitly assessing protein solution rheology (albeit for small proteins), further advancements in computational hardware are still required to permit the application of that type of models outside of supercomputers, as well as to extend their application to biologics of pharmaceutical interest. Finally, the empirical nature and reduced training data sets of most of these computational schemes limits their transferability to different molecules or formulation conditions. Generally, a new set of statistical correlations or CG force field parameters needs to be derived for each new protein system, which can be inefficient and unsuitable for screening over hundreds of drug candidates and/or formulations. Lai et al.⁷⁰ suggested, as a potential solution to this challenge, the use of databases to relate sets of model parameters to viscosity values for linear interpolation.

Protein-protein and protein-excipient interactions

As highlighted in the previous section, protein interactions are the keystone to understanding and predicting the different physical instability phenomena in protein solutions. The balance between different intermolecular forces (from both protein–protein and protein-excipient interactions) and solution conditions defines the likelihood that proteins self-associate, which in turn triggers the different stability issues discussed above. While protein–protein interactions account for the total contribution from attractive (e.g., hydrophobic, van der Waals, dipoles, hydrogen bonding) and repulsive forces (e.g., electrostatic, sterics), the protein environment plays a critical role in constraining the relevant ensemble of configurations that determines the net interactions and protein behavior.¹⁰⁶ As such, at diluted conditions, long-range interactions and the distribution of interacting sites at the protein surface govern net protein interactions due to the large average separation distances between proteins. On the other hand, at elevated concentrations, factors such as multibody interactions, crowding effects, spatial correlations and ion binding become important in controlling the distribution of proteins in solution, and thus they are equally important to net protein interactions.^3,111,200 From an experimental standpoint, protein interactions at diluted conditions are captured through parameters such as the second virial coefficient $B_{22}$ and the dynamic interaction parameter $k_D$, while other parameters like the Kirkwood-Buff integral $G_{22}$ and structure factor $\mathit{S}\left(\mathit{q}\right)$ are used for assessing high-concentration protein interactions.^201,202 Likewise, the preferential interaction parameter ${\mathrm{\Gamma }}_{23}$ is also used as a metric for protein-excipient interactions.²⁰³ In fact, this assortment of parameters is central to many theories, models and heuristics for predicting the high-concentration behavior of protein solutions.^201,204 This section focuses on summarizing the computational approaches used to capture these parameters and their application to evaluate and predict protein stability problems.

Among the different parameters used to describe protein–protein interactions in diluted conditions, $B_{22}$ is arguably the most broadly used. This parameter provides a measurement of the orientational- and solvent-averaged protein–protein interactions, and it has been correlated to phenomena such as protein aggregation¹¹² and protein-phase behavior.¹⁶⁹ One of the earliest, yet widely used computational approaches to calculate $B_{22}$ is derived from the Derjaguin-Landau and Verwey-Overbeek (DLVO) theory for colloidal systems.²⁰⁵ This theory represents proteins as spherical particles interacting through a combination of van der Waals and screened electrostatic forces. An extended version of DLVO (xDLVO) has also been developed to incorporate other types of interactions such as dispersion forces and osmotic potential.^206,207 Different reports have shown that both DLVO and xDLVO can predict the overall trends of $B_{22}$ for protein solutions as a function of ionic excipient concentration²⁰⁸ and in the presence of polymers.²⁰⁶ Nonetheless, these theories provide only qualitative representations for protein interactions, as they are unable to accurately capture how $B_{22}$ is affected by anisotropy in terms of both protein shape and surface heterogeneity.^106,209

Pusara et al.²¹⁰ recently developed a coarse-grained xDLVO model to correct for the spherical approximation on the DLVO theory, though the results only show a modest improvement for predicting $B_{22}$ values as a function of ionic strength on both globular proteins and immunoglobulins. These results are not surprising, as Grünberger et al.²¹¹ previously demonstrated that a CG representation of at least one bead per residue is the minimum resolution-level required to reasonably reproduce steric effects on $B_{22}$ on different classes of proteins. In this regard, the low-resolution CG models from Kim and Hummer⁵⁵ and Blanco et al.¹⁴ offer suitable options to predict $B_{22}$, as they provide a resolution of one CG-site per residue that interacts through amino acid-specific short-range attractions and screened electrostatics (akin to a classical DLVO theory). These models differ by the type of protein interactions they were validated to capture. The former model was designed to reproduce “lock-and-key” and moderately strong protein interactions in order to evaluate the formation and behavior of protein complexes.⁵⁸ On the other hand, the latter model can represent the effects of pH and ionic strength on $B_{22}$, and it has been implemented to evaluate weak protein–protein interactions on globular proteins⁴⁸ and mAb formulations.⁵⁰ Notably, atomistic and high-resolution CG models such as MARTINI have also been implemented for reproducing experimental $B_{22}$ values of globular proteins in different solution conditions, though they require significant re-parameterization of the corresponding force fields.^212,213

Although low-resolution CG models constitute simple but accurate protein representations for evaluating $B_{22}$ and identifying the effect of specific residues in the colloidal stability, they present practical limitations for their implementation in formulation development applications. When considering the use of these models for screening a broad range of formulation conditions, the computational burden of simulating hundreds of conditions makes them unsuitable for this type of application due to typically constrained timelines for drug development. As such, both continuum and simplified CG models provide alternative approaches for efficiently evaluating protein interactions. Wertheim’s perturbation theory has been applied to evaluate $B_{22}$, showing the ability of this type of continuum model to semi-quantitatively capture the relationship between anisotropic protein interactions and the nature of the solution buffer and other ionic excipients.^104,176 Likewise, Roberts and collaborators have extensively evaluated the use of simplified CG models for capturing the behavior of $B_{22}$ on globular proteins⁸⁷ and mAbs.^15,65 In an initial report, Calero-Rubio et al.¹⁵ used a series of mAb model representations (ranging from 1 to 12 CG-sites) to compare the effect of a number of physical factors and model parameters (e.g., hinge flexibility, charge distribution, and the strength of attractions) on protein–protein interactions at diluted and concentrated conditions. The authors found that highly anisotropic charge distributions lead to nonphysically realistic $B_{22}$ values, while hinge flexibility has minimal impact on protein interactions. More importantly, they identified that mAb models of 6 or 12 CG-sites provided a reasonable balance between computational cost and numerical uncertainty when comparing $B_{22}$ results against low-resolution CG models. The same group later fit these simplified CG models against experimental data of $B_{22}$ versus ionic strength for two different mAbs in order to evaluate their ability to predict the behavior of mAb solutions at high-concentrations.^63,64 Interestingly, it was shown that both types of CG models are able to reproduce the osmotic compressibility over a wide range of mAb concentrations as long as net protein interactions are repulsive or mildly attractive. More recently, Shahfar et al.⁶⁵ also compared different simplified CG models against $B_{22}$ values for five different mAbs and arrived at a similar conclusion, where domain-resolution models (e.g., a 6- or 12-bead mAb representation) are only able to reasonably predict $B_{22}$ for repulsive and mildly attractive conditions. The authors also found that a 12-bead mAb model with explicit incorporation of charged amino acids on the protein surface can reproduce $B_{22}$ for conditions dominated by attractive electrostatic interactions.

Computational models are also often used to evaluate the static structure factor $\mathit{S}\left(\mathit{q}\right)$ and the osmotic compressibility as metrics for protein–protein interactions at elevated concentrations.²⁰¹ $\mathit{S}\left(\mathit{q}\right)$ is related to the Fourier Transform of the protein radial distribution function and provides a measurement of spatial correlations at all length-scales.¹⁶ The osmotic compressibility is instead related to the zero-limit value of the structure factor (i.e., $\mathit{S}\left(0\right)$), and therefore it provides information regarding how protein molecules are correlated with each other in the bulk solution.²⁰² From the standpoint of continuum models, there exists a number of different approaches for calculating $\mathit{S}\left(\mathit{q}\right)$ based on integral equation theory, where all of these models represent molecules as spherical particles interacting via various continuous or discontinuous potentials.^82,105,107 As discussed in the previous section, the high-concentration behavior of protein solutions is driven by anisotropic interactions, and thus these isotropic models present limited applicability to capture $\mathit{S}\left(\mathit{q}\right)$ for protein solutions.⁸² Nonetheless, akin to the case of $B_{22}$, Wertheim’s theory has also been extended to evaluate the role of specific anisotropic attractions on $\mathit{S}\left(0\right)$ for globular proteins¹⁰⁴ and mAbs.¹⁹⁹ Alternatively, Minton^214,215 developed a continuum model for evaluating $\mathit{S}\left(0\right)$ (defined by the author as $M_{\mathit{W},\mathit{app}}=M_W\mathit{S}\left(0\right)$) from light scattering experiments, which incorporates molecular anisotropy by representing proteins as hard convex particles based on the approximation than molecular crowding is the dominant force at elevated protein concentrations. Note that this model only considers excluded volume effects as direct protein interactions, while it implicitly captures the effects of longer-range attractions and repulsions by: (1) accounting for a thermodynamic equilibrium between monomer and protein clusters; and (2) allowing the protein diameter to differ from the actual molecular size. As a result, Minton’s model has been successfully applied to evaluate protein self-association at high concentrations,^216,217 and to study how weak protein clusters are related to solution viscosity^201,218,219 and protein-excipient interactions.^220,221 However, due to the inherent assumptions regarding protein morphology and the nature of the intermolecular interactions, this model might lead to overestimation of net protein interactions and cluster formation when the solution behavior is dominated by strong protein-protein/excipient interactions.^201,221,222

Similarly, simplified CG models have been commonly used for fitting experimental $\mathit{S}\left(\mathit{q}\right)$ profiles from small-angle scattering experiments on mAb systems.^69,82,85 For instance, Corbett et al.⁸² evaluated $\mathit{S}\left(\mathit{q}\right)$ for mAb solutions at different pH and protein concentrations as high as 160 mg/mL, using a 3 CG-site mAb model that interacts only via short-range attractions. The authors found that such a simplified CG model is able to capture generic features of $\mathit{S}\left(\mathit{q}\right)$ related to the quaternary protein structure and the nearest-neighbor interaction shell of the proteins. However, when the 3 CG-site mAb model is parameterized using solely data from diluted conditions, it is unable to reproduce high-concentration protein interactions as measured by both $\mathit{S}\left(\mathit{q}\right)$ and $\mathit{S}\left(0\right)$. Likewise, Dear et al.⁸⁵ and Chowdhury et al.^88,198 used a 12 CG-site model to evaluate experimental $\mathit{S}\left(\mathit{q}\right)$ from two different mAbs and a bovine immunoglobulin in various formulation conditions, where the CG-sites interact through a weak short-range attraction and an electrostatic repulsion. Although these reports show that the 12 CG-site model is able to reasonably reproduce the regions of $\mathit{S}\left(\mathit{q}\right)$ related to the bulk behavior and nearest-neighbor shell (i.e., the low- and intermediate-$\mathit{q}$ regions), the addition of a strong attractive potential to the outer CG-sites of the model is required for fitting the data from net attractive formulations.

Most of the above in-silico models mainly focus on evaluating protein–protein interactions, while the effects of excipient are implicitly treated and absorbed in the resulting model parameters. That is, a different set of parameters for the corresponding CG force field might be required for representing different solution conditions of the same protein. As a result, these models efficiently screen protein–protein interactions for different formulation conditions, but they are unable to provide information about the mechanism of action for a given excipient. In that regard, atomistic models have been used to explore how different excipients interact with the protein surface and/or disrupt protein-protein/solvent interactions.^{21,204,223,224} Generally, these studies of protein-excipient interactions follow similar docking methodologies than those used for identifying binding free-energies for the formation of protein-ligand complexes.¹⁰ However, unlike protein–ligand interactions, pharmaceutically relevant excipients such as carbohydrates, nonionic surfactants and free amino acids interact weakly with proteins (Figure 4). As a result, these studies have shown that protein-excipient interactions occur through multiple interacting regions rather than specific binding pockets. Trout group proposed the use of the preferential interaction parameter ${\mathrm{\Gamma }}_{23}$ as a metric for comparing protein-excipient interactions, which quantifies the excess of excipient molecules in the vicinity of the protein as compared to the bulk solution.²⁰⁴ Cloutier et al.^30,143 have investigated ${\mathrm{\Gamma }}_{23}$ for the interaction of three different mAbs with sugars (sorbitol, sucrose and trehalose) and ionic excipients (NaCl, arginine and proline) in the context of protein aggregation and solution viscosity. The authors found that carbohydrates often interact with aromatic residues, whereas free amino acids interact through both charge–charge and cation–π interactions. Nonetheless, the effects of these excipients on protein stability are not generalized, as they lead to different aggregation and viscosity behaviors across the different tested mAbs. The same authors further extended this approach by developing a machine learning algorithm that combines the calculation of ${\mathrm{\Gamma }}_{23}$ with in-silico molecular descriptors (e.g., protein charge, surface area, hydrophobicity) in order to screen biotherapeutic formulations.²²⁵

Figure 4. Illustrative example of protein-excipient interactions for the variable region of a mAb as calculated by the preferential interaction parameter (${\mathrm{\Gamma }}_{23}$). Panels show the interactions of the antibody with different excipients: (a) Proline; (b) arginine-HCl; and (c) NaCl. Coloring indicates local values of ${\mathrm{\Gamma }}_{23}$, where red indicates preferential inclusion (i.e., attractive interactions). Notably, for all excipients, multiple regions of preferential inclusion are identified along the protein surface. Figure adapted from Cloutier et al.¹⁴³

Using a different computational approach, Jo et al.²²⁴ and Somani et al.²²⁶ recently studied the relation between protein-excipient interactions and mAb solution viscosity by evaluating the binding free-energies of a number of excipients (including free amino acids and sugars) to three different Fabs, for which they used the Site-Identification by Ligand Competitive Saturation (SILCS) technology.²²⁷ SILCS consists of simulating the protein in an aqueous environment containing a number of different probe, organic molecules, which provides information regarding the affinity of the protein toward specific functional groups. This information is thereafter used to estimate the binding free-energy of different protein regions to the excipients of interest via a perturbation approach, and thus reducing the overall computational cost of individually simulating each of the protein-excipient systems. By applying SILCS, the authors were able to map the regions of the proteins with higher affinity toward interactions with the different excipients, where some of these regions exclusively interact with specific excipients. Comparison of the patterns for protein-excipient interactions against similar patterns for protein–protein interactions and experimental data of high-concentration viscosity yielded an interesting negative correlation between the number of binding sites for a given excipient and viscosity.²²⁴ The analysis led the authors to hypothesize the mechanism of action by which lysine increased viscosity for one of the studied mAbs.²²⁶ However, the computational results by SILCS could not explain how some of the excipients such as arginine modulate the behavior of the studied mAbs. While SILCS presents a promising methodology for efficiently performing excipient selection during formulation screening, further efforts are still needed to adapt this new technology to the type of problems faced in the biotherapeutic industry, as it was also acknowledged by the same authors.

Future directions

This article reviewed computational protein models of different resolutions, as well as a selection of in-silico studies applied to elucidate some of the biopharmaceutically relevant physical stability issues in high-concentration protein solutions. Models developed to assess protein aggregation,^73,100,144 protein-phase separation,^83,177,178 elevated solution viscosity,^81,195 and protein interactions were included.^85,225 Table 4 highlights several of the computational models used for studying physical instabilities in proteins, which are described throughout the review. Modulating these problems constitutes one of the main goals when developing phase-appropriate drug product formulations in order to gain understanding about their underlying mechanisms and the role that protein environment plays on these phenomena. Besides the inherent challenges in studying these complex issues, experimental assessment of concentrated protein solutions is often restricted by limitations in instrument capabilities, material availability and accelerated timelines for development. In this regard, computational models like those reviewed here could facilitate the identification of suitable drug formulation by helping analyze and augment the information drawn from experiments, providing molecular insight about the stability of the protein in different formulation conditions, and reducing the number of formulation parameters to be assessed in-vitro. Clearly, efforts toward that direction have already begun in earnest, as significant advancements have been achieved in terms of developing novel models and methodologies for efficiently studying the behavior of pharmaceutical proteins such as mAbs, as well as for predicting how protein stability is affected by formulation parameters (e.g., pH, ionic strength, excipients). However, this review also highlights many of the challenges to fully integrate modeling approaches in formulation development workflows.

Table 4. Summary of computational models for studying physical instabilities in protein solutions

Model	Resolution	Underlying Approach	References
Protein Self-Association and Aggregation
Lumry-Eyring Nucleated-Polymerization	Continuum	Aggregation kinetics model that accounts for protein folding, nucleation, and growth by monomer addition and aggregate coalescence Applied to aggregation kinetics data of globular proteins and mAbs Limited applicability at elevated protein concentration	^Citation100,Citation117,Citation118
Arosio et al.	Continuum	Aggregation kinetics model that accounts for protein folding, nucleation, and growth by monomer addition and aggregate coalescence Aggregate size and morphology are considered through Smoluchowski’s population balance for both diffusion- and reaction-limited aggregation Accounts for the effect of aggregate size on viscosity to apply the model at elevated protein concentrations	^Citation109,Citation120,Citation228
Michaels et al.	Continuum	Aggregation kinetics model that includes primary nucleation by oligomerization, secondary oligomerization by aggregate fragmentation, and growth by fibril elongation Applied to amyloidogenic proteins Not applicable to aggregates with morphology beyond fibrils	^Citation101,Citation115,Citation122,Citation123
Baftizadeh et al.	Atomistic	All-atom model for poly-valine and Aβ peptide using AMBER99 force field in Bias-Exchange Metadynamics Studied the conformational transition from disordered oligomers into parallel β-sheets	^Citation145,Citation229
Tofoleanu et al.	Atomistic	All-atom model for different variants of the islet amyloid polypeptide using CHARMM36 force field in NPT-MD simulations Evaluated the dynamics and conformational stability of fibrils	^Citation149
Schwierz et al.	Atomistic	All-atom model for Aβ peptide using CHARMM27 force field in NPT-MD simulations Evaluated the secondary-nucleation aggregation mechanism and the role of solvation during fibril elongation	^Citation151,Citation152
Nichols et al.	Atomistic	All-atom model for different FV variants of a mAb using CHARMM27 force field in NPT-MD simulations Atomistic simulation used for evaluating local structural dynamics, which were combined with calculation of molecular descriptors and APR predictions Applied to interpret experimental outcomes related to solution viscosity and aggregate content	^Citation154
Rojas et al.	Low-resolution CG	UNRES model applied to Aβ peptide using NVT-MD simulations Studied the effects of sequence on the folding, nucleation, and elongation stages during fibrillation	^Citation62,Citation230,Citation231
Vàcha and Frenkel	Simplified CG	Peptides represented as spherocylinders existing in a two-state conformation to distinguish between soluble and aggregation-prone monomer conformations; model evaluated via MC simulations Studied the early stages of the fibrillation mechanism of peptides, including primary and secondary nucleation and growth	^Citation71,Citation73,Citation123
Phase Separation
Kern and Frenkel	Simplified CG	CG model for patchy particles; patches represented as solid angles on the surface Commonly used for representing globular proteins with different potentials for short-range interactions Applied to evaluate the effect of anisotropic interaction on self-association and phase behavior of patchy molecules (including proteins and soft matter)	^Citation83,Citation171–173,Citation232
Kastelic et al.	Continuum	Based on Wertheim’s perturbation theory for representing anisotropic interactions in proteins; patches are modeled as interacting points on the protein surface Applied to globular proteins and mAbs to evaluate the effect of directional interactions on phase behavior, aggregation, viscosity, and protein-excipient interactions	^Citation104,Citation175,Citation176,Citation178–180,Citation199,Citation233
Sun et al.	Simplified CG	Represents mAbs with 3–7 beads and incorporates flexibility at the hinge region; Simulations carried out with Discrete-MD Evaluated the effect of mAb morphology and structure flexibility on the critical point and shape of binodal curve for LLPS	^Citation177
Dignon et al.	Low-resolution CG	Represents intrinsically disordered proteins with 1 bead per residue, incorporating backbone flexibility through harmonic spring for two-body interactions Applied for studying how phase behavior of disordered proteins is affected by structure dynamics, specific interactions, and post-translational modifications	^Citation57,Citation80,Citation164
Transport Properties
Mode-Coupling Theory	Continuum	First-principle theory for the prediction of rheological properties that arise from density fluctuations in colloidal glasses Applied to characterize the concentration-dependent viscosity of globular proteins solutions MCT fits well protein solutions governed by long-range repulsions, but it fails to reproduce the rheological behavior of proteins intermediate-range order due to short-range attractions and long-range repulsions	^Citation118,Citation189,Citation191,Citation193,Citation194
von Bülow et al.	Atomistic	All-atom model for several small globular proteins (no larger than 14 kDa) using AMBER99SB*-ILDN-Q force field in NPT-MD simulations Used for calculating the effect of protein concentration and cluster formation on solution viscosity	^Citation195
Chaudhri et al.	Simplified CG	12- and 26-bead mAb representations; beads interact with each other via site-specific short-range attractions and long-range electrostatics; model used in NPT-MD simulations Studied the solution structure of up to six different mAbs to evaluate the formation of antibody clusters, networks, and higher-order structures, which were correlated with the viscosity behavior; electrostatic attractions between Fabs and the Fc were identified as the driving force for mAb self-association	^Citation66–68
Wang et al.	Simplified CG	12-bead mAb representations, where beads interact with each other via short-range attractions and long-range electrostatics; model used in Brownian Dynamics simulations with an implicit model for hydrodynamic forces Explicit calculations of solution structure, diffusivities and viscosity of several mAbs; the formation of weakly interacting clusters was identified as the cause for changes in transport properties	^Citation69,Citation70
Lai et al.	Atomistic	All-atom mAb model using CHARMM36m force field in NPT-MD simulations; model applied for 27 different mAbs MD simulations were used for relaxing the protein structure and evaluating ensemble-averaged values for several different descriptors such as solvent-accessible surface area, SAP and surface charge map values In-silico results were combined with experimental viscosity data in a Machine Learning feature selection approach	^Citation144
Dear et al.	Simplified CG	12-bead mAb representations, where beads interact with each other via short-range attractions and long-range electrostatics; model was parameterized from fitting experimental SAXS data to calculated $\mathit{S}\left(\mathit{q}\right)$ from NVT-MD simulations Applied to evaluate cluster size distribution at different formulation conditions and protein concentrations for two mAbs and a polyclonal IgG; results from cluster distributions were used to predict concentration-dependent viscosity	^Citation85,Citation88,Citation198
Protein Interactions
Blanco et al.	Low-resolution CG	1-bead per residue representation, where beads interact with each other via residue-specific short-range attractions and long-range electrostatics; model was parameterized from fits to experimental $B_{22}$ values Applied to studying dimerization of globular proteins and assessing the effect of formulation conditions on diluted protein interaction (from $B_{22}$) for antibodies and globular proteins	^Citation14,Citation48–50
Calero-Rubio et al.	Simplified CG	1 to 12 beads per protein representation, where beads interact with each other via residue-specific short-range attractions and long-range electrostatics; model was implemented in Monte Carlo simulations Applied to studying the effect of solution conditions on both diluted and concentrated protein interaction (from $B_{22}$ and $\mathit{S}\left(0\right)$, respectively) Authors identified that a 6-bead mAb representation is the minimum resolution for accurately predict $B_{22}$ at repulsive and mildly attractive conditions For strongly attractive interactions, it was required to explicitly account for surface charged residues in the model to reproduce $B_{22}$	^Citation34,Citation63–65,Citation87
Corbett et al.	Simplified CG	3-bead per mAb representation, where beads interact with each other via short-range attractions Model was used to characterize $\mathit{S}\left(\mathit{q}\right)$ profiles from SAXS experiments for a mAb in different solution conditions and protein concentrations Force field parameters required to be concentration-dependent in order to reproduce both $\mathit{S}\left(\mathit{q}\right)$ and $\mathit{S}\left(0\right)$	^Citation82
Minton	Continuum	Proteins represented as hard convex particles; model explicitly accounts for thermodynamic equilibrium between monomers and protein clusters Developed for characterizing $\mathit{S}\left(0\right)$ from light scattering experiments Applied to evaluate protein self-association at high concentrations and the relation between protein clusters and solution viscosity and protein-excipient interactions	^Citation201,Citation214,Citation215,Citation217,Citation218,Citation220
Shukla et al.	Atomistic	All-atom protein model for assessing protein-excipient interactions using various versions of the CHARMM force field in NPT-MD simulations ${\mathrm{\Gamma }}_{23}$ is calculated from simulations based on the distribution of both solvent and excipient molecules with respect to the protein Applied to characterize the interaction surface of globular proteins and mAbs with different cosolutes and pharmaceutically relevant excipients A machine learning algorithm was recently developed to combine the calculation of ${\mathrm{\Gamma }}_{23}$ with in-silico molecular descriptors (e.g., protein charge, surface area, hydrophobicity) for screening excipients against aggregation propensity and rheological behavior of mAb solutions	^Citation30,Citation143,Citation204,Citation225

The most pressing limitation is without a doubt the lack of benchmarking of existing computational models against diverse experimental data sets, which should contain information about solution behavior over multiple proteins for a wide range of formulation parameters. Many of the current protein models are borrowed or adapted from other fields (e.g., amyloidosis, crystallization, colloids), where similar stability issues are also of interest. However, such models and their corresponding force fields are unable to accurately capture the behavior of proteins in pharamaceutically relevant conditions (e.g., non-biological pH or buffers, interactions of proteins with surfactants). Likewise, for models with a high level of atomistic detail, explicit simulation of concentrated protein solution remains computationally prohibitive, and thus they generally rely on statistical correlations for linking the simulation outcomes to experimental behavior. Due to the lack of comprehensive data, the conclusions of those models present limited transferability and extrapolability to proteins/formulations outside the range of experimental data used to build such correlations. This latter point is also by itself another challenge.

Notably, many in-silico methodologies highlighted here have led to a seemingly evident conclusion regarding the unique nature of protein behavior with respect to formulation conditions. That is, even for proteins sharing a high level of homology between each other (e.g., mAbs of the same isotype), the effects of an excipient, for instance, might drastically vary from one protein to another. As such, there is a need to develop and implement advanced statistical or mathematical methodologies that can evaluate complex data sets and contextualize the broad, but apparently contradictory effects of a given formulation across different proteins. In this regard, Trout’s group has already taken some steps in this direction by combining, via machine learning algorithms, the results from different atomistic models for predicting various instability behaviors with experimental data.^144,225 Lastly, even with the use of simplified protein models, the computational time for simulating concentrated protein systems remains unfeasible for efficiently screening over hundreds of formulation conditions and/or drug candidates. Consequently, further efforts are needed regarding the development or implementation of complementary theories or numerical approximations that facilitate reducing either the computational time or the number of simulations. Methodologies such as those proposed by Mahynski et al.¹⁸⁵ and Hatch et al.²³⁴ for extrapolating modeling results for predicting phase separation and $B_{22}$, or the theories developed by Carmichael and Shell²³⁵ and Jin et al.¹⁸³ for multi-scaling in aggregation and phase separation problems might provide a solution to this latter challenge. The complexity of the phenomena related to protein stability and high-concentration formulations is likely to continue fueling the advancement of many more computational models and methodologies, facilitating the continued significance of in-silico workflows to the development of biologic drug products.

Abbreviations

Table

APR	Aggregation-prone region
BSA	Bovine serum albumin
CDR	Complementarity-determining region
CFD	Computational fluid dynamics
CG	Coarse-grain
DLVO	Derjaguin-Landau and Verwey-Overbeek theory
$k_D$	Dynamic interaction parameter
xDLVO	Extended version of DLVO
Fab	Fragment antigen-binding region
FC	Fragment crystallizable region
FV	Variable fragment
IgG	Immunoglobulin G
$G_{22}$	Kirkwood-Buff integral
LENP	Lumry-Eyring Nucleated-Polymerization model
LLPS	Liquid-liquid phase separation
mAb	Monoclonal antibody
MCT	Mode-coupling theory
MD	Molecular Dynamics
MW	Molecular weight
MC	Monte Carlo
$B_{22}$	Osmotic second viral coefficient
${\mathrm{\Gamma }}_{23}$	Preferential interaction parameter
SILCS	Site-identification by ligand competitive saturation
SAXS	Small-angle X-ray scattering
$\mathit{S}\left(\mathit{q}\right)$	Static structure factor
NVT	Thermodynamic canonical ensemble
NPT	Thermodynamic isothermal-isobaric ensemble

Acknowledgments

The author thanks numerous MSD colleagues for their input at various stages of this work. Dr. Hanmi Xi and Dr. Timothy Rhodes are gratefully thanked for comments to the manuscript. Dr. Pouria Mistani is also thanked for many helpful and stimulating discussions.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

The author(s) reported there is no funding associated with the work featured in this article.