2.1. Representative Interactions of Cellular Crowding

The current picture of the inside of cells is that of a crowded environment, as the one shown in Figure ​Figure11a, where biochemical reactions have to take place in a highly concentrated solution with various species of finite sizes, which includes both low molecular weight (molecular crowding) and high molecular weight (macromolecular crowding). Note that the term “crowding agents” is commonly used for molecules similar in size to the considered protein. The influence of water molecules, which are much smaller, is considered as solvation,⁴⁰ while for much bigger molecules (such as macromolecules that can be considered immobile obstacles), the effect is referred to as confinement.⁴¹ This can be understood by picturing the cell as a room full of people: while the motion of a person is influenced by bystanders, smaller entities (like ants or dust) do not change it. The presence of voluminous furniture would delimit the room but not be part of the flow. Historically the effects of crowding agents have been divided into soft chemical/physical interactions and volume exclusion (or hard-core repulsion).⁴²⁻⁴⁵

Open in a separate window

Figure 1

Schematic diagram of crowding inside a cell and of its effect across protein binding. (a) Scientific illustration of the macromolecules inside a cell of Escherichia coli, inspired by the figure presented by Goodsell³⁹ with the cytoplasm in blue and purple and the cell membrane in yellow. The magnified portion is a magnification of the cytoplasm constituents. (b) Outline of the effects of crowding on protein folding, diffusion, and binding (from top to bottom) inside the cell. Crowding effects can be divided into volume exclusion and soft interactions. The former has a repulsive nature, which tends to enhance the stability of folding and binding and decrease the translational diffusion, as summarized in the first column. Soft interactions can instead be classified as attractive or repulsive (second and third column, respectively). Attractive soft interactions tend to counteract volume exclusion, which has a destabilizing effect on both folding and binding. Conversely, repulsive soft interactions have a stabilizing effect. Both interactions hinder rotational and translational diffusion.

The term “excluded volume” was used by Minton, who coined the expression “macromolecular crowding” in 1981,² to indicate the effect that volume exclusion has on the energetics and transport properties of molecules in a highly volume-occupied medium.⁴⁶ Volume exclusion is a steric effect arising from the impenetrable nature of atoms that can not overlap because of the hard-core or van der Waals repulsion and only involves molecules rearrangement and, thus, affects the entropic component of protein stability. In this setting, the activity of a molecule, for example, a protein in equilibrium between the unfolded state and native folded state, is raised when crowding molecules are added because its excluded volume is inaccessible to the centers of the other crowding molecules.^46,47 This favors its folded state because it corresponds to a smaller excluded volume. For the same reason, the oligomerization of the protein is promoted, and the association (dissociation) constant for ligand binding is increased (decreased). In general, according to this view, crowding facilitates processes leading to a reduction in excluded volume, like protein folding, oligomerization, complexation, aggregation, and condensation.⁴⁸

To characterize the effects of crowding given by this definition, the past decades have seen an explosion of experiments using large artificial polymers, such as Ficoll, PEG, and dextran to recreate a volume exclusion effect and allow the generalization of effects that may depend on the specific protein crowder. Even if Ficoll and PEG are still often considered as noninteracting with proteins,^49,50 many studies have argued that Ficoll and PEG show weak chemical interactions with the test protein.⁵¹⁻⁵⁵ In a recent computational study on the behavior of PEG and Ficoll as crowders, Ostrowska et al.⁵³ investigated how they influence a specific enzyme called NS3/4A. By employing atomistic simulations, they found that while the enzyme’s affinity for its substrate remains similar with or without the crowding agents, the speed at which it performs its catalytic function decreases with PEG and increases with Ficoll. This suggests that these crowding agents may affect the enzyme’s function through specific interactions rather than the more generic volume exclusion effect. The study also showed that both crowding agents made contact with the enzyme and slowed down its movement. Additionally, they found that the crowding agents influenced the structure of a protein component called NS4A, which caused it to adopt helical structures. Overall, PEG had a slightly stronger interaction with NS3/4A, while Ficoll formed more hydrogen bonds with the enzyme. This result was in line with previous experimental measures.⁵⁵

Although for more than 30 years crowding theories have emphasized steric repulsion, in the last decades, the number of studies probing the competition between the excluded volume effects and soft interactions has been continuously rising.^43,56−62 Soft chemical/physical interactions change the enthalpy of proteins and include water–protein and solute–protein interactions. In this context, some of the key weak interactions are (i) van der Waals interactions, (ii) hydrogen bonds, (iii) ionic interactions, (iv) dipole–dipole interactions, and (v) weak noncovalent hydrophobic interactions. Most soft interactions are attractive and are expected to favor expanded conformations (e.g., the unfolded state of a protein will expose sites for attractive interactions, such as hydrogen bonding and hydrophobic contacts), thereby leading to destabilization.⁴⁵ Thus, they are typically believed to oppose and reverse the hard-core repulsion effect of crowding. However, among soft interactions, there is a strong repulsive one arising from the opposition of charges with the same sign. This repulsion has a stabilizing impact similar to that of the hard-core effect:⁶³ in recent years, it has been experimentally and theoretically shown that chemical interactions can be both stabilizing or destabilizing.⁶⁴⁻⁶⁶ The acknowledgment of soft interactions in theories of crowding has brought the conclusion that inert polymers—used in most experiments or simulations aimed at assessing hard-core repulsion effects—are not good mimics of the in-cell environment.^67,68 To distinguish between these inert polymers and the biomolecules that can, indeed, be found in the cell, we will use the terms “inert crowders” and “crowders,” respectively. To investigate the crowding effects in cell-like conditions, different types of “crowders” with their chemical–physical characteristics have to be considered. The importance of considering more complex systems to obtain a realistic representation of the cellular environment and investigating both the entropic and enthalpic effects that take place in vivo has been recently reviewed by Pastore et al.⁶⁹

2.2. Protein Interactions in the Crowded Cellular Environment

Protein–protein interactions are essential for protein function and cellular pathways formation. Since these interactions are also involved in the development of diseases, they are important targets for drug design and the artificial design of protein complexes. Most studies usually distinguish between non functional and functional partners to focus on the latter. “Functional partners” can be defined as proteins that have been evolutionarily selected to have functional interactions with a few partners⁷⁰ and avoid interactions with thousands of other types of macromolecules in the cell. However, crowding makes “nonfunctional partners” interact despite this evolutionary selection, which leads to a competition between functional and nonfunctional binders.^{26,38,70−72} This competition implies that to predict functional interactions, weak and potentially nonfunctional ones have to be characterized, as well.^26,73

In addition to this, in the context of the cellular interior, complexity of the classification, itself, has been challenged.^42,74−76 For example, sequential metabolic enzymes seem to have evolved to associate weakly,⁷⁷ and similarly, weakly bound functional protein complexes are observed in metabolic, regulatory, and signaling pathways. Sukenik et al.⁷⁸ hypothesized that such complexes create an additional fuzzy network of weak interactions in the cell, which gives the cell the ability to detect, signal, and/or directly initiate regulatory processes quickly and effectively in response to external stresses or internal signals (such as rapid volume variation). According to them, the cell environment is finely tuned to optimize weak interactions networks among its protein machinery. In a case like this, it is difficult to see the involved proteins as functional partners or crowders that are weakly interacting with each other; in a crowded environment, all interactions have to be considered.

The interaction between two molecules is defined as a correlation between their positions, orientations, and motions.⁷⁹ This can happen even when the two partners are distant because of the presence of electrostatic or solvent-mediated interactions. These long-range interactions can be nonspecific or specific.⁷⁹ Nonspecific interactions do not lead to the formation of a specific complex and are, thus, more likely to be independent of the partners’ orientation. Specific interactions are strong and lead the diffusional search of binding partners for one another followed by colocalization and recognition of the compatible binding sites. This determines the formation of a complex with a defined structure and interface. Because of the specific structure that has to be achieved, interactions are expected to depend on the mutual orientations of the partners. In both cases, the importance of a match between different physicochemical properties at a global level and not only between the interacting sites (if any) has been shown: for the partners approach, their total charge, isoelectric point, hydration free energy, and total electrostatic energy have to match.⁸⁰⁻⁸² Significant commonalities for the pH of maximal stability⁸³ were found, as well. Simplified model potential of mean forces, including these general features of the overall intermolecular interaction, have been shown to obtain semiquantitative agreement with experimental measurements of the osmotic second virial coefficient⁸⁴ and of concentration-dependent light scattering and osmotic pressure.⁷⁹

These results underline the importance of taking into account the crowded cellular environment in which the investigated proteins live to understand how their binding can be influenced by these nonspecific interactions.

3.1. Predicting if Two Proteins Interact

Many existing techniques for predicting protein–protein interactions and their networks have been developed using only sequence information^{117−119,151,157−167} since they can be applied to proteins for which structural information is unknown or that are intrinsically disordered. Most sequence-based approaches start from the hypothesis that pairs of proteins that are similar to pairs of interacting proteins have a higher chance to interact. Different algorithms have been implemented to identify similar regions.

For instance, Shen et al.¹⁶⁸ developed a learning algorithm based on a support vector machine, along with a kernel function and a conjoint triad feature for amino acid description, and trained the model using over 16 000 diverse pairs. In this respect, the effectiveness of current sequence-based methods in constructing interaction networks and understanding disease mechanisms is highly dependent on the availability and reliability of training data (as reviewed by Murakami et al.¹⁶⁹), which means that predictive methods, especially those based on parameter training (i.e., supervised machine learning methods), need a large amount of data to improve their predictive ability.

The importance of the data set construction has been addressed by Li et al.,¹¹³ who proposed another sequence-based prediction method called Scoring PRotein INTeractions (SPRINT). This algorithm identifies possible binding partners regions that are similar to sites on known complexes with a multiple spaced-seed approach and then eliminates elements that occur too often to be involved in interactions. SPRINT was shown to be several orders of magnitude faster than other state-of-the-art sequence-based programs, including PIPE2¹¹⁶ and two machine learning (ML)-based methods proposed by Ding et al.¹⁷⁰ and Martin et al.¹⁷¹

Another key step for prediction methods based on ML is the representation of the input data. Romero-Molina et al.¹¹⁸ developed a procedure that transforms pairs of amino acid sequences into a ML-friendly vector whose elements represent numerical descriptors of residues in proteins. This numerical encoding method was then implemented by the same group to develop a support vector machine model, PPI-Detect.¹¹⁸

In the same year, Yao et al.¹¹⁹ proposed a new residue representation method named Res2vec and combined it with deep learning techniques to predict interactions starting only from sequence information (DeepFE-PPI).

Other methods based on sequence information for predicting residues involved in intermolecular binding were developed using coevolution.¹¹⁵ The idea is that many proteins have evolved to form specific molecular complexes, and the specificity of their interactions plays a crucial role in their proper functioning. Therefore, these interactions impose constraints on the protein sequences because the network of inter-residue contacts must be maintained. It is plausible to assume that any sequence change accumulated during the evolution of one interacting protein is counterbalanced by changes in the other protein.

Pazos et al. developed a method for detecting correlated changes in multiple sequence alignments, which is useful to a collection of interacting protein domains.¹⁷² The results revealed that positions exhibiting correlated changes in both interacting molecules tend to be near the protein–protein interfaces. This observation opened up the possibility of applying statistical inference techniques based on the maximum entropy principle^173−178 to predict pairs of residues that come into contact solely based on their sequence information^179−181 and, more specifically, in their evolutionary process.

Many methods for interaction prediction focus on protein–protein interaction networks since genome-scale data for different species have been rapidly increasing. By starting from the network of an organism and evaluating the network topology features of proteins, these methods can predict new interactions.¹⁸² In the past, protein–protein interaction networks were thought of as static, defined by interactions constant in any cell location and condition.¹⁸³ The current picture is that interactions are dynamic and that their occurrence depends on a set of conditions, such as spatial or temporal variations.¹⁸⁴ In other words, these networks are now modeled as dynamic systems. Computational approaches can come in aid of experiments in the prediction of the variation and dynamics of these networks in the context of crowding.²⁶In silico methods tend to detect protein complexes by associating them with meaningful clusters in the protein–protein interactions network. They can be classified into three categories:¹⁸⁵ cluster-quality-based,^186,187 node-affinity-based,^188,189 and ensemble clustering methods.¹⁹⁰ The first one uses clusters as measuring units, while node-affinity-based methods measure the affinity between nodes inside clusters. Finally, ensemble clustering methods combine the clusters selected by different methods to mine final complexes.

However, all these methods rely on experimentally determined protein complexes data sets, which are afflicted by three main problems: (i) they are incomplete,¹⁹¹ (ii) they contain many false positives caused by experimental conditions,¹⁹² and (iii) they have been measured in conditions that do not replicate the cellular environment.

Other techniques look at the functional similarity between proteins using Gene Ontology (GO) terms.¹⁹³ GO terms give information about a protein’s localization within the cells, participation in biological processes, and associated molecular functions.¹⁹³ Interacting proteins belong to the same pathway and, thus, tend to participate in similar processes and/or have similar functions and/or live in similar cellular compartments.^194−196 Therefore, these methods can classify interactions with a similarity score of GO terms.^196−198

3.2. Predicting Binding Interfaces and Poses

Binding takes place on portions of the protein interfaces (or binding sites); for natural dimeric interfaces, they range from a change in solvent-accessible surface area of 850–10 000 Ų (up to 7000 Ų for heterodimers).¹⁹⁹ These protein regions have been vastly characterized, for example, for what concerns their amino acid composition.^200,201

Studies of protein–protein interactions can be divided into two categories: (i) those identifying only the interfaces and (ii) those also predicting the complexes 3D structures.

Many methods in the former category only require the sequence of the investigated proteins because structural studies have shown that binding regions are characterized by a combination of geometrical and chemical complementarities^{200,202−204} that can be predicted from the proteins’ amino acid sequences.^205,206 In this context, homology modeling is very effective, as well. An example is BSpred,¹²⁰ which is based on neural networks. The algorithm underwent extensive training using sequence-based features, such as protein sequence profile, secondary structure prediction, and hydrophobicity scales of amino acids. Other recent sequence-based methods for predicting binding sites include simplified long short-term memory (LSTM) network,²⁰⁷ the DeepPPISP method that uses a combination of local and global sequence features,¹²⁴ and the DELPHI method, which combines CNN and recurrent neural network (RNN) in an ensemble structure.¹²⁶

Other protocols take as input the proteins’ structures and focus on features that are known to characterize interacting molecular surfaces: the preferentially hydrophobic composition of the binding interfaces or the role of van der Waals interactions.^{107,200,208,209} Milanetti et al.^210,211 proposed the use of 2D Zernike polynomials to rapidly compare protein molecular iso-electron density surfaces and assess their shape complementarity. The protocol called Zepyros (ZErnike Polynomials analYsis of pROtein Shapes) has been developed and distributed freely as a web application by Miotto et al.¹²⁵ and applied to characterize different biological systems^212−214 and to optimize biomolecule interfaces.^215−218 The work by Grassmann et al.²¹⁹ extended the procedure to account for electrostatic similarity.

Gainza et al.¹²⁸ proposed that proteins involved in similar interactions might possess shared fingerprints, regardless of their evolutionary origins. They introduced MaSIF (molecular surface interaction fingerprinting), a conceptual framework utilizing geometric deep-learning techniques to capture crucial fingerprints relevant to specific biomolecular interactions, including protein pocket–ligand interactions and protein–protein interaction. Two years later, the same framework was used to design novel protein interactions.²²⁰

Several other methods for predicting protein binding sites and residue–residue contacts involve a mixed approach between sequence information and structural information and are based on ML approaches. PAIRpred, a Support Vector Machine approach developed by Asfar Minhas et al.,¹²¹ employs a specific pairwise kernel and a combination of structural and sequence-based features, such as position-specific score matrices (PSSMs), position-specific frequency matrices (PSFMs), and solvent accessibility predictions. The same features used in PAIRpred are used in a graph convolutional neural network (GCNN) proposed by Fout et al.¹²²

Sequence-based and structural features are used also in BIPSPI (xgBoost Interface Prediction of Specific-Partner Interactions), which was developed by Sanchez-Garcia et al.¹²³ More recently, Krapp et al.¹²⁹ introduced PeSTO, which acts directly on protein atoms without the need for parametrization of the system’s physics and is provided as a web server.

Techniques based on the structure of the investigated proteins have become powerful tools for predicting protein–protein association because of the success of structure modeling methods; when protein structures are not available from experiments, it is now possible to predict them with a certain degree of reliability starting from the sequence information. The recent advancements in the fields of protein structure and complex predictions have been recently reviewed by Wodak et al.²²¹ The state of the art in structure modeling is assessed by the Critical Assessment of Structure Prediction (CASP).²²² The top-ranked methods in the latest competitions of CASP have reached excellent results. Some software for structure prediction are Rosetta²²³ and D-I-TASSER.²²⁴ Rosetta developers introduced an all-atom force field that focuses on short-range interactions (e.g., van der Waals, hydrogen bonds, and desolvation effects) and discards long-range electrostatics. Predicted structures are refined via a Metropolis Monte Carlo-based refinement protocol. This method reaches a Cα root-mean-square deviation (RMSD) = ~1.5 Å.²²³ A more recent method is D-I-TASSER, which can predict protein structure and function and was extended from the I-TASSER method developed by the Zhang lab,²²⁵ which integrates threading and deep learning. Its pipeline was ranked in the first place in the CASP15 experiment in all categories of protein structure prediction.²²⁶

Still, the most popular and effective tool for protein structure prediction is AlphaFold, the top-ranked method in the 13th CASP edition in 2018.²²⁷ Two years later, in CASP14, AlphaFold2 was presented²²⁸ and again outclassed the other methods. AlphaFold2 was able to predict the structure of a protein only from its amino acid sequence (we shall remind the reader that the human genome has been mapped since 2001²²⁹) with a median global distance test (GDT) score of 92.4 across all targets.²³⁰ The efficiency of this tool in predicting the coordinates of the heavy atoms of the structure relies on the integration of novel neural network architectures and training procedures based on the evolutionary, physical, and geometric constraints of protein structures. In terms of RMSD₉₅ (i.e., the Cα root-mean-square-deviation at 95% residue coverage), it allowed a median backbone prediction accuracy of 0.96 Å against the 2.8 Å achieved by the next best-ranked prediction method. As for all-atom predictions, its accuracy was 1.5 Å RMSD₉₅ compared with the 3.5 Å RMSD₉₅ of the next best-performing approach.²²⁸

Prediction of the interfaces is not always sufficient; to understand the physical mechanisms involving protein complexes, determination of their 3D structures is a critical step. AlphaFold2 has given impressive results even in this context: when applied as a method for complex pose prediction, it exceeded docking approaches for heterodimers prediction,^139,231 even if the modeling of antibody–antigen complexes was unsuccessful (11% success²³¹). Other limitations were still not resolved: AlphaFold2 can not predict intrinsically disordered proteins or regions.²³² Moreover, it is based on a multiple sequence alignment and, thus, fails in the prediction of novel structures. In October 2021, Alphafold2 was extended to multiple chains and took the name of AlphaFold-Multimer.²³³ AlphaFold-Multimer was able to predict homodimers even better than heterodimers, despite still failing in the prediction of antigen–antibodies binding. AlphaFold-Multimer can predict heteromeric interfaces with DockQ²³⁴ higher than 0.23 and 0.8 in 70% and 26% of the cases, respectively.²³³ Compared with the performance of AlphaFold, the improvement was +27 and +14 percentage points in each case.²³³ Homomeric interfaces were predicted in 72% of the cases with a high accuracy 36% of the time (+8 and +7 percentage points, respectively).²³³ This year, Liu et al.²³⁵ proposed a further enhancement of AlphaFold2-Multimer based on a novel quaternary structure prediction system, MULTICOM.²³⁶ MULTICOM improves AlphaFold-Multimer predictions by generating diverse multiple sequence alignments and structural templates using both sequence and structure alignments and combining the AlphaFold-Multimer confidence score with the complementary pairwise model similarity score to rank models. It then classifies them through multiple complementary metrics and refines the models via a Foldseek²³⁷ structure-alignment-based refinement.

Many improvements were implemented on AlphaFold-Multimers, which resulted in a better performance compared with existing approaches, but its accuracy in complex prediction was still lower than AlphaFold2’s accuracy for tertiary structure prediction.²³⁵

Despite the limitation given by the need to have experimentally resolved homologous structures, methods like AlphaFold that predict protein–protein interactions using sequence homology are among the most promising approaches.²³⁸ Given the importance of the structural characterization of protein complexes, it is not surprising that the number of large complexes deposited each year in the Protein Data Bank (PDB) is growing rapidly.²³⁹ A significant contribution to this trend comes from the advancement of experimental technologies, such as the developments of methods based on electron microscopy (EM), as reviewed by Saibil et al.,²⁴⁰ which are particularly suitable for systems of macromolecular assemblies. Consequently, new computational approaches based on homology information have been developed.^241−245

Baspinar et al. proposed PRISM,¹³⁴ which combines structural similarity and accounts for evolutionary conservation in the template interfaces. One of the most used homology-based methods is SWISS-MODEL, which is a pioneered, automated protein homology modeling server that has undergone consistent enhancements over the past 20 years.^{138,246−249} Recently, its modeling capabilities have been expanded to incorporate the modeling of both homo- and heteromeric complexes by utilizing the amino acid sequences of the interacting partners as the initial reference.¹³⁸ Nevertheless, these template-based approaches can offer a satisfactory prediction only when homolog complexes are available.

If the structures of the investigated molecules are available, useful computational tools are docking algorithms, that can take in input two (or more) component proteins and assemble them into putative models of the protein complex. As thoroughly reviewed by Pagadala et al.,²⁵⁰ in the last two decades, more than 60 docking tools have been proposed. In the following, we will provide a quick overview of the most efficient ones. The performances of docking algorithms are usually tested in competitions like CAPRI.²⁵¹

The docking methods participating in CAPRI can be classified into two approaches: the first one uses a matching technique that describes the partner as complementary surfaces,¹³³ and the second one simulates the actual docking process and computes the pairwise interaction energies.^252−254 The former predicts the static structure of protein complexes, while the latter is based on simulations that also explore the kinetics leading to the binding and will thus be discussed in Section 3.3.

In complementarity-based docking, efficient search strategies are fundamental for global docking software to achieve good performance because of the high cost of the search stage. Many available global docking software rely on FFT correlation search algorithms that enable full systematic searches through translational and rotational degrees of freedom. The FFT approach was originally introduced to evaluate only the shape complementarity²⁵⁵ but has been later expanded to electrostatic interactions (e.g., FTDock²⁰²) or to both electrostatic and solvent terms, for example, in ZDOCK.²⁵⁶ ZDOCK, which also considers atomic contact energy to estimate electrostatic corrections and further improve its results, is one of the best-performing FFT-based docking programs.¹³³ However, simpler scoring functions can reduce the computational costs: MEGADOCK¹³⁷ is a FFT-based rigid-body docking tool similar to ZDOCK²⁵⁶ but ~7.5 times faster thanks to the evaluation of only shape complementarity and electrostatic for the scoring.²⁵⁷ Better results have also been obtained with PIPER,¹³¹ the first FFT-based approach to use pairwise structure-based potentials (extracted from structures of protein–protein complexes) called Decoys As the Reference State (DARS).²⁵⁸ Its protocol starts with the generation of a large decoy set of docked conformations that are used as a reference state. Using DARS, the contact frequency between two specific atom types in the native state and in the decoys are compared. By rewarding the occurrence in the interface of the atom pairs that frequently interact in the native complexes, PIPER produced at least 50% more hits in 19 of the 33 tests for enzyme–inhibitor complexes compared with ZDOCK.¹³¹ While still obtaining more hits than ZDOCK in 12 of the 16 test problems, PIPER was less accurate for antigen–antibody.¹³¹ This limitation was addressed six years after its development with the introduction of a new nonsymmetric potential extracted from antibody–protein complexes.²⁵⁹

All the presented docking methods only accept structures as input. Yan et al. proposed HDOCK,¹³⁵ which can accept both sequence and structure inputs. HDOCK is a hybrid docking strategy that combines a FFT-based global docking with template-based modeling. By looking at the available homologous complexes in the protein data bank, it can incorporate binding interface information -if any- into traditional global docking.

Among rigid body docking algorithms, there are also a few global programs that employ alternative search strategies, such as randomized search or local shape matching, across the entire protein structure, like PatchDock¹³⁰ and SymmDock.¹³⁰ The former uses geometric hashing algorithms to perform a global protein–protein docking using local shape descriptors (i.e., surface patches) to predict protein–protein and protein–small molecule complexes. The latter can only predict the structure of a homomultimer with cyclic symmetry given the structure of the monomeric unit.

A common problem of all these docking algorithms is that they can overlook plausible complex structures because of how the poses are ranked: after the search, for a pair of input proteins, the typical docking program generates tens of thousands of different binding poses, including both near-native (i.e., almost correct) and incorrect ones. Good scoring functions are fundamental for identifying the best poses among the proposed set. Even if many types of scoring functions have been developed, they often simplify the interactions on the basis of the shape complementarity and treat solvation effects implicitly.²⁶⁰ Even more advanced programs, like PatchDock, which utilizes geometric hashing for rapid surface patch matching, and ZDOCK, which also employs an advanced pairwise method to efficiently incorporate shape complementarity, do not consider many physical components, including entropy contributions from solvent molecules and proteins. In addition to this, Wass et al.²⁶¹ showed that for interacting proteins even the incorrect poses have relatively favorable scores, although they are lower than the correct ones.

Kozakov et al.¹³⁶ developed the ClusPro docking server, which improves the docking prediction by taking into account entropic effects.²⁶² The protocol starts by implementing PIPER as a rigid body docking step and then performs a RMSD clustering of the 1000 lowest energy structures generated. To choose the most likely near-native clusters, their size is evaluated instead of their energy values: the centers of the largest clusters are selected rather than the lowest energy structures. The selected structures are then refined using energy minimization. A similar procedure is followed by MEGADOCK,¹³⁷ which is based on a FFT-based rigid-docking scheme. MEGADOCK classifies protein pairs as interacting if among the docking-generated poses there are clusters of similar poses with significantly favorable docking scores compared with the others.

All the presented methods allow for only partial protein flexibility. Methods that perform side-chain optimization are able to reach higher-accuracy solutions compared with more rigid methods.²⁶² Desta et al.²⁶² recently tested some of the most used docking servers on the same data set. They showed that ZDOCK and pyDock²⁶³ (another rigid body server similar to ClusPro) obtained the worst results compared with SwarmDock,²⁶⁴ a flexible docking algorithm that was able to produce slightly more models in the top 1 prediction (predicting good models for 23.5% of the 51 targets). This method models flexibility using a linear combination of elastic network normal modes and combines a local search with a variation of the particle swarm optimization (PSO) algorithm.²⁶⁵ In this algorithm, each particle represents a potential solution and moves in a D-dimensional space (including Cartesian space, orientational space as quaternions, and coefficient space for a linear combination of anisotropic elastic network normal modes). During the search, each particle adjusts its position according to its own experience and the swarm’s experience. PSO can simultaneously optimize translational, orientational, and conformational degrees of freedom. It is worth mentioning that SwarmDock was implemented in 2010 in a study on macromolecular crowding⁸⁷ that will be discussed in the next sections. However, the SwarmDock server¹³² has the limitation of handling only files with less than 10 000 atoms.

Another flexible docking protocol, HADDOCK2, was proposed by Van Zundert et al.⁹¹ and has now become one of the most exploited data-driven approaches. Data-driven methods use biochemical and/or biophysical interaction data (e.g., chemical shift perturbation data resulting from NMR titration experiments²⁶⁶) to guide the docking process. In HADDOCK2, experimental data are used to define the simulated annealing in torsion angle space to introduce partial flexibility. Thanks to this flexibility, HADDOCK2 has been shown to yield more accurate models than ClusPro.²⁶²

However, considering only side-chain flexibility has been argued to be insufficient: in the CAPRI experiment, the majority of failures were attributed to the inaccurate prediction of conformational changes in proteins during protein–protein interactions.²⁶⁷

3.3. Predicting the Structural Determinants of Binding

Different theories to describe the kinetics leading to the binding pose have been proposed. The prevailing view has evolved from the early “lock-and-key” hypothesis to the “induced-fit” and “conformational selection” models.²⁶⁸ In the former, partners undergo very small changes of backbone and side-chain atoms upon binding (root-mean-square deviations of 0.6 and 1.7 Å, respectively²⁶⁹). However, it is known that proteins are dynamic: both “induced-fit and “conformational selection” hypotheses take this aspect into account. According to the former, the interactions between two approaching structures induce conformational changes. Conversely, the latter states that the protein’s bound state can be explored even in the absence of the molecular partner.

Potential transition states have been deeply investigated by experimental studies, like double mutant cycles and paramagnetic relaxation enhancement.^270−272 However, the resulting data are often indirect or limited to specific cases (like metalloproteins or proteins with attached paramagnetic spin labels). Computational methods can instead provide atomic-level observations of the association pathways and explored conformation for a diverse set of protein complexes.

The most employed methods to accomplish this purpose are all-atom (AA) molecular dynamics (MD) simulations, which provide both structural and dynamical insight into protein–protein binding. Some of the most widely used simulation software packages include CHARMM,¹⁴² AMBER,¹⁴¹ GROMOS,²⁷³ and GROMACS.¹⁴⁰

This approach allows us to witness the evolution of the system made of the involved protein partners in terms of association/dissociation events along the simulated trajectory. In this way, it is possible to sample a large number of binding poses while gaining information about intermediate states characterizing the association pathway. Thus, MD simulations are often at the core of simulation-based docking techniques.^274,275

In this context, two intermolecular interaction potentials describing Coulombic and van der Waals interactions are present between atoms of the two molecular partners to modulate their interaction. In particular, given two atoms a_l and a_m holding partial charges q_l and q_m, the Coulombic interaction between them can be computed as

1

where r_lm is the distance between the two atoms, and ϵ₀ is the vacuum permittivity. van der Waals interactions can instead be calculated by, for example, the Lennard-Jones potential:

2

where ϵ_l and ϵ_m are the depths of the potential wells of a_l and a_m, respectively, and R_min^l and R_min^m are the distances at which the potentials reach their minima.

Commonly available computational power for MD simulation typically allows them to sample only a few binding poses in the short microsecond simulation frame given the long time scale of the interconversion between poses.

Longer time scales (as well as larger systems) can be achieved with coarse-grained (CG) models²⁷⁶ in which several atoms are lumped together in effective interaction sites (beads). Compared with AA representations, CG models have fewer interactions, smoother free energy profiles, and fewer degrees of freedom. Even if CG representations usually need additional terms to preserve the native structure that limits their conformational sampling, alternative structural restraints have been proposed.²⁷⁷

To find a compromise between the benefits and limitations of AA and CG models, many multiple-resolution modeling strategies have been developed. Some of them treat at atomistic resolution the biomolecule’s functionally relevant portion,²⁷⁸ while the rest of the system is described with a CG model, which results in a lower computational cost compared with AA model.²⁷⁹ One of the most used force fields is MARTINI.^280,281 As a downsize, these methods often require lengthy reference all-atom simulations and/or the usage of off-shelf coarse-grained force fields to parametrize the CG-modeled part.

The treatment of the interactions between different levels of resolutions has to be determined, as well.¹⁴⁷ To overcome these limitations, Fiorentini et al.¹⁴⁷ proposed a novel multiresolution modeling scheme for proteins, dubbed coarse-grained anisotropic network model for variable resolution simulations (CANVAS). CANVAS, which can be implemented on the most commonly used MD platforms, allows for a user-defined modulation of the resolution level throughout the system structure and a fast parametrization of the potential without reference simulations.

In combination with CG and AA models, structure-based models (SBMs), also called G-type models, can be used.^282,283 While CG methods model the coarse-grained resolution of the system, SBMs measure the coarse-grained resolution of its force field.

They usually take an experimentally solved native conformation of the protein to determine the native contacts and use them to define the interaction potential for the simulation by introducing bias toward the native state. This reduces the force field complexity without loss of essential information according to energy landscape theory and the principle of minimal frustration,²⁸⁴ which states that a protein energy landscape has a funnel-like shape biased toward the native state. This ensures the dominance of native over non-native interactions, thereby enabling an efficient folding. The prevalence of native interactions in the folding event has been confirmed by atomistic simulations.²⁸⁵ Thanks to this simplified interaction potential and the consequent reduced computational costs, structure-based models are often used to simulate many cycles of folding and unfolding.

Another popular approach—that can use both atomic and coarse-grained representations—to shed light on molecular processes at time scales that are not accessible with single unbiased MD simulations is the combination of MD simulation with Markov-state modeling (MSM). Popular software packages are pyEMMA¹⁴³ or MSMBuilder.¹⁴⁶

Another approach is speeding up the sampling of rare events by enhancing the conformational sampling performed during the simulation, usually by scaling the temperature or using collective variables (CVs). CVs are identified by dimensionality reduction methods that can determine which properties are of interest for the description of the high-dimensional conformational ensembles.²⁸⁶ In addition to facilitating the interpretation of the huge amount of data produced during a simulation, CVs can direct the conformational sampling to efficiently cover the free energy landscape.²⁸⁷

Some of the algorithms using CVs to enhance the sampling include steered MD,²⁸⁸ umbrella sampling,²⁸⁹ Hamiltonian replica exchange MD (REMD),²⁹⁰ and metadynamics.²⁹¹

Alternatively, the methods using temperature to enhance the sampling include simulated annealing,²⁹² simulated tempering,²⁹³ temperature replica exchange MD (T-REMD),²⁹⁴ and replica exchange with solute tempering (REST²⁹⁵ and REST2²⁹⁶).

In all these cases, a central role is given to the force field, whose quality determines whether the sampled conformations are realistic. Because of this, force fields are constantly being improved.

The most popular atomic-level protein force fields family include AMBER,²⁹⁷ CHARMM,²⁹⁸ OPLS,²⁹⁹ and GROMOS.³⁰⁰ They explicitly model each atom or, in the case of GROMOS, all heavy and nonaliphatic hydrogen atoms, and describe the interactions between bonded and nonbonded atoms with mathematical functions whose parameters are fitted to structural or thermodynamic quantum mechanical or experimental data.

Among the AMBER force fields, ff99SB³⁰¹ has been often recommended for proteins. It has been shown to improve amino-acid-dependent properties and reproduce the differences in amino-acid-specific Ramachandram Map using amino-acid-specific CMAPS.³⁰² Used with the OPC (optimal point charges) water model,³⁰³ it was also found to improve the accuracy of atomistic simulations of disordered proteins.³⁰⁴ However, with time, ff99SB showed weaknesses in some amino acid side chain dihedral parameters;³⁰⁵ in 2019, an updated model with improved backbone profiles for all 20 amino acids, named ff19SB,³⁰⁶ was proposed.

Similar backbone CMAP potentials were used two years prior for CHARMMM36m,²⁹⁸ an updated version of CHARMM22,³⁰⁷ which had been for a long time the first choice among CHARMM force fields for proteins simulations. This improvement increased the accuracy of generating ensembles for disordered proteins.

The evaluation of force fields treatment of disordered proteins is commonly done by comparing them^308,309 because the intrinsic complexity and “roughness” of the disordered proteins energy landscape can reveal force field deficiencies and strengths. For instance, a comparison between CHARMM36m and ff19SB concerning their ability to describe disordered proteins while retaining proper description of the folded ones has been recently performed by Abriata et al.³¹⁰ Both CHARMM36m and ff19SB were applied with the recommended water model, respectively modified by TIP3P²⁹⁸ and OPC.³⁰⁶ CHARMM36m was found to favor more compact states than ff19SB-OPC, which leads to the overstabilization of aggregates and secondary structures. ff19SB better predicted weak dimerization of a soluble protein while still being able to reproduce the expected aggregation of another system. Conversely, CHARMM36m predicted residue-wise alpha helical propensities better than ff19SB.

In the last decades of protein simulations, other popular force fields have been OPLS/AA³¹¹ and its modification OPLS-AA/L.³¹² In 2015, Robertson et al.³¹³ developed OPLS-AA/M, which proposes a new parametrization of the backbone and side chain dihedral. With a similar strategy, Harder et al.²⁹⁹ presented OPLS3 a year later. The former has an improved ability to simulate disordered proteins, while the latter well simulates protein–ligand binding.

A similar reparameterization against experimental data was proposed for the GROMOS force field.³⁰⁰ Even if it resulted in an improved prediction of secondary structure propensities and structural parameters, investigations on how these new parameters could perform for intrinsically disordered proteins are still ongoing.

4.1. Effect of Crowding on the Single Protein Structure

Protein dynamics link protein structure and function: to meet and interact with other proteins and perform their function, proteins move and explore their conformational space.^405−407 Thanks to in vitro(408) and in silico(409) experiments and theoretical models,⁴¹⁰ it has been known for decades that protein folding free energy landscapes are shallow; free energy differences between folded and denatured states and activation free energies are relatively small (0–40 kJ/mol).³⁶⁶ This means that the stability of a protein can be influenced by variations in environmental conditions. Protein folding and stability have been well studied under simplified conditions, mainly through in vitro experiments, but it is yet unclear how these protein properties can be modulated by the crowded interior of live cells.⁴¹¹

The effects of crowding on conformational variation have been often described as a shift in equilibrium between native and unfolded states, but the observation of a non-native state that was different from unfolding suggested a more complicated scenario^{66,359,365,366,393,412} in which entropic and enthalpic effects may (co)exist. Harada et al.⁸⁵ hypothesized that, in crowded conditions, entropic constraints should explain the disfavoring of (thermally or chemically) unfolded states, while enthalpic contributions could promote alternative non-native states.^{66,359,365,366,393}

The general take-home message of in silico studies that only consider the entropic excluded volume effect is that it enhances stabilization of the native state by increasing the free energy of the unfolded configurations for both single- and multidomain proteins.³⁹⁵ The stabilization of the more compact native states by volume exclusion is consistent with many experiments.^49,396,413

Some have measured an enthalpy-driven stabilization, as well; for example, Benton et al.⁴¹⁴ affirmed that different-sized crowders have the same effect on a protein. This stabilization has been presumed to derive from specific interactions with nearby crowder protein surfaces,^365,366,415 such as edge-to-edge β-sheet interactions,^366,416 or altered solvent properties.⁴¹⁷ More generically, enthalpic effects have been connected to a destabilization of the folded forms given by protein–protein nonspecific interactions. These interactions counteract the stabilization of the folded forms, which results in a reversed overall effect of the crowding,^418,419 even if sometimes a stabilizing total effect has been observed both in vitro(415) and in silico.^63,365 The importance of soft interactions has been discussed in many in silico studies that had to add them to the volume exclusion to match experimental values.^71,420 These contrasting results could depend on the fact that the crowding effect on protein folding has been shown to depend on a variety of factors, including crowders’ occupancy, shape, and size,^50,98,394 as well as the considered protein.^98,421

The inconsistency found in the literature is, of course, also determined by the many different models of crowding that can be chosen. A simple example is given by the evolution of the description of crowding effects on the 20-residue-long Trp-cage protein shown in Figure ​Figure88. Trp-cage, whose 3D structure was obtained with NMR in 2002,⁴⁰² is a valuable model system for understanding protein folding dynamics: because of its reduced length and fast folding,⁴⁰³ it has been often tested both inside and outside the crowding field. Figure ​Figure88 shows on a timeline some of the studies that have investigated Trp-cage folding in the presence (on top) and in the absence (on bottom) of crowders.^422−424 One year after the insertion of its 3D structure in the PDB,⁴²⁵ Zhou⁴⁰³ explored the folding free energy landscape of Trp-cage in explicit solvent but without crowders. He observed a two-step folding mechanism with a metastable intermediate state that justifies protein’s fast folding. At 300 K, 70% of the population corresponded to the folded state. One of the first investigations of the effect of crowders on Trp-cage was performed in 2010 by Tsao et al.³⁶⁵ They introduced an implementation of replica-exchange simulations with discrete MD (in which the simulations proceed according to ballistic equations of motion) to study the effect of spherical crowders on an atomistic description of Trp-cage. They only investigated the volume exclusion effect and found that the protein is pushed to adopt a compact state regardless of whether it is the native conformation or not. Two years later, Predeus et al.³⁵⁹ increased the realism of the description by trying to achieve a balance between realistic description and computational efficiency. They introduced a three-component multiscale modeling scheme in which proteins were modeled at an atomistic level, the crowder proteins were modeled at a coarse-grained level (using the PRIMO coarse-grained model⁴²⁶), and the surrounding aqueous solvent was represented as an implicit solvent (using a GB-type implicit solvent model). By combining this model with temperature replica exchange sampling, the simulations enabled the exploration of relatively long time scales that reached into the microsecond range to investigate the impact of protein G crowders on Trp-cage. This approach allowed them to go from simple volume exclusion considerations to including charges and LJ terms. While Tsao et al., by using spherical crowders, indicated a destabilization of extended and unfolded conformations, the interactions between crowders and Trp-cage led to significant populations of partially folded structures that were not observed in corresponding simulations without crowders. This suggests that crowding determines not only a stability shift for Trp-cage but also essential conformational changes. With a dielectric constant of ϵ = 80, Predeus et al. observed two near-native conformations in agreement with what was found by Tsao et al.³⁶⁵ But while Tsao et al. found that the radius of gyration of the second state was larger than that of the native state, for Predeus et al., it was smaller.

Open in a separate window

Figure 8

Schematic representation of six different computational studies of the protein Trp-cage structure and dynamics. On top, three articles (from left to right: (359, 365, and 393)) that show an example of how the computational investigation of the dynamic in crowded conditions of a case study protein, in this case Trp-cage, has evolved over time. For each paper, the most relevant bibliographic data are reported together with a representative illustration, and the year of publication is marked on the timeline. The same information is shown for the papers on the bottom. The study on the left refers to the insertion of the Trp-cage structure in the PDB databank;⁴⁰² the other two (from left to right: (403 and 404)) show the evolution of computational methods that investigated its dynamic without considering crowding.

A similar system was simulated in 2015 by Bille et al.:³⁹³ they chose a system size close to the one studied by Predeus et al. by using the B1 domain of protein G (GB1) and bovine pancreatic trypsin inhibitors (BPTI), which have similar size (but different net charge) compared with protein G, as crowders. Here, crowders were represented in atomistic details, even if with the assumption of a restricted internal dynamic limiting them to the folded state. In this study, they employed Monte Carlo (MC) methods and an AA representation with an implicit solvent force field. As in the previous study, they found that crowders hinder the protein from adopting its global native fold, but this time the interactions between Trp-cage and the crowding agent were found to be specific and involving a select few key residues, therefore denying the generality of the nature of Trp-cage–crowder interactions.

The discrepancy between the findings of Predeus et al. and Bille et al. may, indeed, be partially attributed to the different net charges of BPTI (+6) and GB1 (−4). Since Trp-cage carries a positive net charge (+1), there are more potential attractive electrostatic residue pair interactions between Trp-cage and GB1 than between Trp-cage and BPTI, which likely contributes to the less specific Trp-cage–crowder interactions in the GB1 case. Additionally, while two prolines in BPTI were identified as key residues for interaction with Trp-cage, GB1 lacks any prolines. Therefore, there are at least two reasons to anticipate differences in the nature of Trp-cage–BPTI and Trp-cage–GB1 interactions.

The comparison between these two crowding studies highlights the importance of selecting appropriate crowders to accurately represent the crowded cellular environment.

Trp-cage has also been used as a model system to test the Replica exchange with hybrid tempering (REHT) method proposed by Appadurai et al.⁴⁰⁴ REHT is a parallel tempering technique to allow efficient and accurate conformational sampling by accelerating water dynamics. Compared with the standard REST2²⁹⁶ method, it guided Trp-cage to the experimentally measured native structure in ~100 ns instead of 300 ns.

Even if in this study crowders were not considered, the improvement of the mapping of the free energy landscapes is an important requirement for simulating crowded conditions. Indeed, when dealing with crowded environments, the exploration of the proteins conformational landscape is slowed by their significant size and by the conformational and dynamical variation of both crowders and proteins. To tackle this challenge, one can either apply this kind of technical improvement to a crowded system or simulate a simplified environment. The latter solution was chosen by the methods presented in the next section.

4.1.1. At the Beginning of Time There Was a Sphere

The complexity of the detailed protein–crowders interactions is unanimously recognized, but many consider volume exclusion to be the primary effect of crowding. In this context, understanding the excluded volume effects is a necessary step to interpret both experiments and all-atom simulations and to develop more complex models. In general, when considering the stability of single proteins or complexes, one may quantify the effect of crowders by computing the difference in free energy (G) of the process in the presence and absence of the crowded environment, i.e., ΔΔG = Δμ_F – Δμ_I, where F and I refer to a generic final and initial state (e.g., complex and monomers, or folded and unfolded protein). We consider a system in which the only varying thermodynamics quantity is the number of molecules so that differences in G are linked to differences in the chemical potential.

To make analytical progress in this framework, scaled particle theory (SPT)⁴²⁷ is often implemented. Worked out by Reiss, Frisch, and Lebowitz in 1959,⁴²⁸ SPT is an equilibrium theory of hard-sphere fluids giving an approximate expression for the equation of state of hard-sphere mixtures and their thermodynamic properties.⁴²⁹ In SPT, the cellular environment is described as a solution filled with crowders and proteins represented as hard spheres. Molecules are considered to be interacting only when touching, and the interacting energy has only a repulsive component that is strongly dependent on distance. The theory provides an estimate of the free energy cost of creating a rigid cavity in the solution on the basis of the ratio (hence the term “scaled”) of the gyration radius of the molecule of interest corresponding to the size of the cavity and the radius of gyration of the crowding agents.⁴³⁰

An approximate analytical solution for Δμ is given by

7

where a and V are the radius and volume of the protein, S and S_c are the surface area of the protein and the crowders, and C and ϕ are the number density and volume fraction of the crowders.³⁶⁹ According to this model, macromolecular crowding exerts big effects on the folding free energy, thereby predicting a monotonic increase in stability as the crowder packing increases.

Despite assuming sphericity for both protein and crowders, this theory was often qualitatively consistent with results coming from molecular simulations of more complex systems composed of spherical crowders and CG proteins.^50,431 Since the conclusions of both SPT theory and simulations of inert spherical crowders were often confirmed by experimental measurements (with inert crowders),⁵⁰ spherical crowders have gained popularity as a swift and straightforward approach and have found widespread application in various types of simulations.^{63,86,100,365,394,395} Because of the limited number of atoms in these systems, the simulations become comparatively rapid, particularly when employing an implicit water model.

However, this approach does not consider chemical interactions, approximates the treatment of solvation, and has an extreme sensitivity to sphere size and density.

4.2. Kinetics of the Binding and Reaction Rates

Diffusion and kinetics of diffusion-controlled reactions are fundamental transport mechanisms that play a crucial role in describing a wide range of chemical and biological processes.⁴⁵⁴ In crowded environments, the presence of other molecules and structures can affect transport properties and, consequently, the dynamics of chemical reactions in different ways. Various techniques, such as “single-particle tracking,”⁴⁵⁵ “fluorescence recovery after photobleaching (FRAP),”⁴⁵⁶ and “fluorescence correlation spectroscopy,”⁴⁵⁷ have been used to measure the diffusion constants of macromolecules in the cytoplasm, as well as in membranes. All experiments show that the in vivo diffusion of proteins is greatly reduced compared with dilute conditions. For in-depth information about the findings of each of these experimental methods, we refer the reader to other reviews by Prindle et al.,⁴⁵⁸ Cai et al.,⁴⁵⁹ and Yu et al.,⁴⁶⁰ respectively. Many computational studies have been performed to capture the complex interplay between crowding, diffusion, and reaction rates: it is by now well established that crowding can significantly affect the diffusion coefficients of macromolecules,^1,461,462 the diffusion-controlled reaction rates,⁴⁶³ and cause shifts in chemical equilibria.⁶³ At membrane interfaces, it has been observed with computer simulations⁴⁶⁴ that diffusion is enhanced because of protein depletion. However, diffusion rates provided by simulations under crowded conditions^10,24,390 have shown that in the cytoplasm diffusion is slowed down by up to a factor of 10. This evidence has also been confirmed by experimental measurements.^465,466 While the most pronounced effects are observed for macromolecules, crowding impacts the dynamics of small molecules, as well. A widely accepted reason for this is the volume excluded by the crowders that limits free diffusion. However, it has been discussed how the excluded volume effects, alone, could not account for the factor of ~10 reduction.^467,468

A factor whose influence on diffusion has been extensively discussed is hydrodynamic interactions (see, for instance, the review by Zegarra et al.⁴⁶⁹), that is, the interactions between molecules or particles mediated by the fluid in which they are immersed. In crowded environments, the presence of other molecules and structures can significantly modify these interactions. CG simulations and hard-sphere models have often addressed their effect in a regime of short-time diffusion (on nanosecond time scales) where hydrodynamic interactions dominate over negligible protein collisions.⁴⁷⁰ Even if some studies have hypothesized that hydrodynamic interactions play only a minor role^373,374,389 in short-time diffusion,³⁶¹ other works have argued that they have a central role⁴⁶⁵ and that they even slow down long-time diffusion rates.^38,470,471 To better investigate their effect, various methods have been developed to incorporate hydrodynamic interactions into MD simulations, for example, Stokesian dynamics. These methods are computationally expensive but highly accurate and are able to model the interactions between solute and solvent molecules. Careful validation of the simulations through comparison with experimental data and/or other numerical methods²⁴ remains important.

The crowding influence of diffusion on a longer time scale has been investigated by recent works based on large-scale atomistic simulations that were able to analyze protein–protein interactions, and it was concluded that the determinant cause for diffusion hindrance is random nonspecific diffusive encounters.^{89,372,389,390} These interactions, called “sticking,” are typically weak and short-lived but have an important effect because of their sheer numbers.⁴⁷² Both in-cell NMR⁴⁷² and computational models³⁷² of the E. coli cytoplasm have explored this topic. Sticking can result in the formation of multiprotein complexes through nonspecific contacts lasting less than 1 μs. The formation of protein clusters may be the cause for the retarded rotational diffusion compared with translational diffusion,^373,374 something that could not otherwise be explained using hard-sphere models. The influence of protein sticking grows for longer time transport dynamics;³⁷² by achieving 5 μs of simulation, Rickard et al.³⁷² measured a smaller diffusion coefficient than other shorter all-atom studies.^373,389

4.2.1. Rolling in the Cell

When only considering the volume exclusion effect of crowding, diffusion is expected to be slowed down. Different experiments at different volume occupation fractions were found to be described by the phenomenological expression:

8

where κ is a free fitting parameter that was found in diffusion experiments in the presence of crowders for the protein carbon monoxide hemoglobin to be κ = 0.36.⁴⁷³ A similar expression can be applied to the rate k_diff at which the reactants get into close contact, which in turn determines the association rate k_–1. k_–1 can, indeed, be expressed as k_–1 = k + k_diff, where k is the intrinsic association rate once the reactants are in contact. The volume exclusion effect of crowding leads to an association rate given by k_diff() = k_diff( = 0)(1 – )^κ. k_diff( = 0) is the classical Smoluchowski equation describing the association rate without crowders, k_diff( = 0) = 4πDa, where a is the size of the target molecule.⁴⁷⁴ Thus, as the system’s excluded volume increases, k_diff decreases since the diffusion time the reactants need to first encounter becomes large. That said, the overall volume exclusion effect on reaction rates is more complex because once the proteins are in proximity, binding equilibria toward the bound state enhances the overall reaction rate.

One of the primary algorithms to simulate reaction and diffusion in continuous space is the Green’s function reaction dynamics (GFRD)^475,476 represented in Figure ​Figure99a. This exact algorithm models the behavior of a diffusive system by cyclically solving the Green’s functions that propagate the probability distribution. GFRD breaks down the many-body reaction–diffusion problem into one- and two-body problems, which can be analytically solved using Green’s functions.^475,476 These Green’s functions are then utilized to create an event-driven algorithm that allows for large jumps in time and space when particles are far apart. When two particles undergo reaction and diffusion within the same system, higher complexity can be employed to describe their interactions. From simple hard spheres, we can move to two spheres with patches where an analytical solution of the Green’s function is still feasible.⁴⁷⁷ Additionally, within the reaction volume where the reactant molecules are present, it becomes possible to locally describe protein binding as a more detailed process. This enables the construction of a multiscale algorithm that computes the spatial dynamics of proteins as hard spheres over large distances while employing more complex methods to describe binding when the proteins share a reaction volume.⁴⁷⁸

Open in a separate window

Figure 9

Levels of representation employed to study diffusion and reaction rates in a crowded environment. (a) The GFRD algorithm defines a sphere around proteins (orange and green), crowders (violet), and solvent (blue) and computes for each molecule the probability distribution of the time and the position of leaving this volume. Each molecule is thus associated with a next event type and a next event time, which are put and executed in chronological order. Following execution, the molecule is propagated, and a new sphere with new events is determined. (b) In the simplest lattice model, space is discretized into voxels (gray graph), and proteins (orange and green) and crowders (violet) can only move to a neighboring voxel. The solvent is either modeled implicitly or not considered. (c) A more realistic representation can be obtained when all the system components are left free to move, and the solvent is explicitly modeled.

Another possible approach is lattice-based methods, represented in Figure ​Figure99b, where each site can hold one molecule, and crowders are represented as already occupied lattice points. A lattice model was used by Meinecke⁴⁷⁹ to simulate the diffusion dynamics of spherical particles surrounded by inert and inactive crowders of any size and shape. Meinecke found that shape and size considerably determine how strongly diffusion is hampered: small crowders hinder diffusion more than bigger ones. The effect is also stronger for elongated crowders and larger diffusing molecules.

Despite their higher computational cost—GFRD can, for example, be up to 6 orders of magnitude faster—⁴⁸⁰Brownian dynamics (BD) are another commonly used method to study how a crowded environment affects association rates. As outlined in Długosz et al.,⁴⁸¹ they can directly account for crowded conditions, molecule sizes, and shapes, and compute the rates at which molecules come together. Figure ​Figure99c shows a general model in which both proteins and crowders (and solvent molecules, if considered) are free to move.

A stochastic simulation was used by Ridgway et al.⁴⁶⁷ to incorporate the physical dimensions of particles into the diffusion and reaction rates. By combining this model with a proteomic-scale evaluation of protein abundance, researchers could approximate the population and diffusion characteristics of the E. coli cytoplasm. This enabled the study of the volumetric impact of macromolecular crowding on biomolecular diffusion and diffusion-limited reactions. However, the authors observed that this effect could not fully account for retarded protein mobility in vivo, which suggested that other biophysical factors have to be considered.

To further investigate the relation of diffusion with the solvent, environment, hydrodynamic, and attractive interactions have to be taken into account. Other methods based on BD have been employed with this aim since these kinds of simulations allow the consideration of interactions and other features of the cellular environment, as well.⁴⁸¹

The specific effect of volume exclusion and hydrodynamic interactions were investigated by Grimaldo et al.⁴⁷⁰ by simulating a spherical model of an Escherichia coli cytoplasm. The model comprised 15 different macromolecule types at physiological concentrations with each represented by a different-sized sphere. In this way, they took into account the polydispersity of a cellular environment and overcame standard colloidal modeling. This polydispersity was shown to slow down larger macromolecules more strongly than smaller particles via hydrodynamic interactions already at nanosecond time scales before protein collision. The results were confirmed by experimental measurements in cell lysate. The authors also suggested that the slowdown may also affect molecular mobilities and escape rates of long-time processes.

The same E. coli model, studied with BD simulations and the Rotne–Prager–Yamakawa (RPY) tensor, was used by Ando et al.³⁸ to investigate the effect of attractive interactions on diffusion. At first, they only considered a soft repulsive potential and hydrodynamic interactions: without adjustable parameters, the in vivo experimental diffusion constant is reproduced. Thus, they hypothesized that excluded volume effects and hydrodynamic interactions are sufficient to explain the large reduction in the diffusion of macromolecules observed in vivo. Next, they added nonspecific attractive interactions: reduction in diffusivity becomes very sensitive to the macromolecular radius, with the motion of the largest macromolecules being dramatically slowed down. Moreover, long-lived clusters involving the largest macromolecules form if attractions dominate, whereas hydrodynamic interactions give rise to significant, size-independent intermolecular dynamic correlations.

The importance of considering soft interactions was shown by Blanco et al.,⁴⁸² as well. They developed the chain entanglement softened potential (CESP) to study streptavidin diffusion among dextrans of different sizes and concentrations. The CESP model is based on an interparticle interaction potential that includes both hard-core and soft interactions and takes into account the hydrodynamic interaction effects of protein complexes. The experimental long time diffusion coefficient predicted with this potential was shown to be in better agreement with the experimental measured one than the coefficient predicted when only considering hard-core repulsion.

The fact that the experimentally observed slowdown of diffusion can be achieved only when considering attractive interactions in addition to volume exclusion and repulsion was again confirmed by Skora et al.³⁹⁷ In 2020, they considered a mixture of Ficoll70 and double-stranded DNAs (dsDNAs) to investigate the effect of the crowder’s shape/type. Using BD simulations, they observed a stronger slowdown of diffusion by cylindrical double-stranded dsDNAs compared with spherical Ficoll70, which was confirmed by in vitro measurements. For a 60% diffusion reduction, a volume fraction of 5% or 35% was needed, respectively, for dsDNA and Ficoll70. The variation in reduction between the two crowders was not only attributed to the larger volume excluded by the former but also to the differences in the attractive interactions.

4.2.3. Looking Closer: All-Atoms Models

Studying diffusion in a crowded environment simulated at atomistic resolution poses a significant computational challenge. Although the dynamics and underlying mathematics are straightforward, the large number of elements can result in an extremely high computational cost. Moreover, particles can react and bind to membranes or microtubules, thereby changing their effective diffusion constant dynamically over time. Transient interactions between proteins can also result in the formation of clusters, which have been proposed as the dominant crowding effect for reduced rotation diffusion.^373,374 This hypothesis focuses on direct protein–protein contacts, but whether diffusion could be affected to a similar degree by proteins that are nearby even if not in contact is less clear.

To evaluate how close proteins have to come before their diffusion properties are impacted significantly, Nawrocki et al.³⁷⁴ performed MD simulations for over 1 μs for a system of 19 chicken villin headpieces with explicit solvent. They showed that the diffusion slowdown is primarily correlated with direct protein–protein contacts (on a 1–100 ns time scale^373,374) because of the formation of clusters rather than indirect interactions via shared hydration layers.³⁷⁴ Consistently, with the Stokes–Einstein relations for translational (D_t) and rotational (D_r) diffusion,

9

and

10

where k_B is the Boltzmann constant, T is the temperature, r is the particle radius, and η is the viscosity of the surrounding medium, they found that rotational diffusion showed a stronger slowdown.

In the same year, von Bülow et al.³⁸⁹ confirmed that nonspecific protein binding (with a complex lifetime of 1–50 ns) and the formation of clusters account for the high viscosity and slow diffusivity in crowded conditions. They performed 500 ns long atomistic MD simulations of ubiquitin (UBQ), the third IgG-binding domain of protein G (GB3), lysozyme (LYZ), and villin headpiece (VIL) and found that at 200 mg/mL, translational and rotational diffusion decreased by a factor of ~4 and ~6, respectively. They computed the effective hydrodynamic radii using the Stokes–Einstein relations and found that it was consistent with the diffusion variation, thereby proving that this relation remains valid in concentrated protein solutions. Finally, they showed that clustering interactions contributed to ~40% and ~50% of the total slowdown in translational and rotational diffusion, respectively.

The dominating effect of cluster formation on diffusion was also found by Bashardanesh et al.,⁸⁹ who simulated systems of a few species of proteins at an increasing concentration (up to 70% of the mass) with AA MD lasting 1 μs. Even if they found a stronger slowdown than von Bülow et al. (probably due to a higher concentration of protein used in their simulations), they agreed in observing that this effect was particularly noticeable for rotational diffusion. The more noticeable hampering of rotational diffusion was explained by the fact that direct interactions from clustering hinder rotational motion and were related to charge distribution on the surface or weak dispersion interactions. Bashardanesh et al. also noticed that rotational diffusion depends on the protein size, unlike translational diffusion and viscosity: at high concentrations, small proteins were slowed down more than big ones. However, the authors did not propose an explanation.

All the mentioned studies explored time scales of ~0.5–1 μs, while a longer simulation of 5 μs, which reaches a time scale relevant for sampling sticking processes, was performed by Rickard et al.³⁷² Their systems more closely replicated the heterogeneity of the cellular environment, as well. They mimicked a small volume of the E. coli cytoplasm with 2 × 10⁵ atoms, which corresponds to proteins and RNA selected among the most abundant cytoplasmic components. This was the first simulation of a bacterial cytoplasm model that reached atomic precision and a μs time scale. The formation of clusters and the consequent diffusion slowdown were observed here as well, even if the observed diffusion coefficients (~1.4–2 nm²/μs) are smaller than those observed in experiment⁴⁸⁵ and in other atomistic simulations.^373,389 The authors hypothesized that this could be explained by the much shorter time scales explored by those studies, which highlights the potential importance of protein sticking and cluster formation for long-time transport dynamics. Rickard et al. showed that protein–protein contacts are highly dynamic and can be broken and reformed repeatedly with a time scale of contacts on the order of ~1 μs that indeed exceeds the simulation time of previous works.^89,389,464

4.3. Dynamics of Complex formation

After navigating the crowded cellular environment, two (or more) functional partners can bind together. It has been suggested that the binding interfaces of proteins have coevolved with the cytoplasm to avoid nonspecific interactions;^72,372 thus, the influence of crowding would be expected to be negligible at this stage of functional complex formation.

That said, some studies that only considered volume exclusion or repulsive interactions found a modest stabilizing effect.^86,281,369 In this approximation, crowding also seems to favor specific over nonspecific complexes since the former are more tightly bound and, thus, have smaller excluded volumes.⁸⁶ The fact that crowding increases the specificity of protein binding in a crowded cellular environment means that it can have functional effects apart from the relative stability of bound and unbound complexes.

However, adding attractive interactions has resulted in different observations.²⁸¹In vitro experiments⁴⁸⁶ have observed that protein crowders can destabilize protein–peptide complexes because of weak nonspecific interactions between test proteins and crowders. The sticking determined by transient weak interactions can both hinder, via nonfunctional competitive binding,⁵⁴ or promote complex formation.⁸⁷ These interactions can result in transient “quinary structures.” The concept of quinary structures introduces an additional level of protein organization beyond the primary, secondary, tertiary, and quaternary levels. These structures are characterized by functional interactions that exhibit weaker and more transient properties compared with quaternary protein structures. They involve the formation of transient encounter complexes, which serve as the initial step in the assembly of low-affinity and transient protein complexes. These complexes are involved in various biological processes, such as electron transfer, the formation of enzyme–substrate complexes, and the self-association of proteins with weak affinities. The understanding of quinary structures provides valuable insights into the dynamic and diverse nature of protein interactions and their functional implications.⁴⁸⁷

Crowding also has an indirect influence on the dynamics of complex formation by shifting the population of local conformations,³⁹⁰ as discussed in Section 4.1.

4.3.1. Interacting Proteins as Spheres

The models that only take into account the volume exclusion effect between spheres predict a shift of the binding equilibrium of proteins toward the bound state: once the associated state has formed, more volume is available for the crowders, and the total entropy of the associated system increases, which makes the bound state more stable.^30,488 To provide a simple way of depicting the effect of volume exclusion on binding equilibria, a reaction of the type A + B AB, where A and B are reactants, can be considered. The latter reaction is fully reversible and can be characterized by a dissociation constant , where k_off is the unbinding rate, k_on is the reaction rate, and the subscript 0 indicates an absence of crowders in the solution. More generally, the reaction dissociation constant can be defined as K_d,, where is the volume occupation fraction of the crowding agents. Thus, the effect on the dissociation constant is purely thermodynamic and can be characterized by the statistics of the states in which the receptor is either free or occupied. Consider a three-dimensional lattice with Ω lattice sites and a total volume V = Ωl³, where l is the lattice spacing. Set a static target or receptor R, and with L proteins and C crowders diffusing into the lattice with a given diffusion constant, if one partner finds the other, the complex RL is formed. At equilibrium, the probability for the complex to be formed can be calculated as

11

where k_B is the Boltzmann constant, T is the temperature, S_b is the number of states in which the complex is formed, S_ub is the number of states in which the protein is unbound, and E_b is the binding energy. After obtaining S_b and S_ub, P_b can be expressed as

12

One can then show that the change in the dissociation constant as a function of is given by K_d, = K_d,0(1 – )^r, where r is a factor that depends on the crowder size.⁴⁸⁹ Thus, volume exclusion shifts binding equilibria toward the bound state. This shift is expected to result in protein aggregation.⁴⁹⁰

To more deeply investigate the thermodynamic effects of crowding on macromolecular associations, hard-sphere theory can be used in conjunction with Brownian dynamics simulations. This has been done by Ando et al.,³⁹⁴ who were the first to simulate a cell-like inhomogeneous system instead of assuming a uniform crowder size, even if the dependence of crowding effects on the composition of the crowders had already been investigated in the past.⁴⁴¹ They mimicked the physiological concentrations of macromolecules in the cytoplasm of Mycoplasma genitalium by using 41 different types of macromolecules represented by spheres with different radii. They found, both with theoretical calculation and simulations, that the heterogeneity of sizes decreases the effect of crowding when compared with the use of uniform crowders with the same mean size at a fixed volume fraction of the system. They also confirmed that smaller crowders show the strongest effects. However, these kinds of simplified studies are limited by the exclusion of shapes, chemical properties, interactions, and flexibility of macromolecules. Ando et al. recognized that the stabilization of the native or complex states of macromolecules predicted by their model is in contrast with what was observed by other experimental and computational studies.^{43,85,414,418}

The fact that considering only volume exclusion results in a wrong prediction of complex stabilization has been recently investigated by Pradhan et al.²⁸¹ They simulated the dimerization of model system GB1 in the presence of lysozyme crowders at two different resolutions by assuming both protein and crowders as spherical beads and then retaining residue-specific structural details. With the former, they found stabilization of the dimers in the presence of the lysozyme crowder, in accordance with the SPT model. The latter residue-specific model showed that soft interactions destabilized GB1 dimerization and were in agreement with the free energy variation obtained from experimental observation. The computed value of free energy variation, despite having the same sign, was much larger than the experimental one. This was hypothesized to be due to the overestimation of the protein–protein interaction strength in the employed Martini 3.0 force field.⁴⁹¹

This research was partially funded by grants from ERC-2019-Synergy Grant (ASTRA, no. 855923), EIC-2022-PathfinderOpen (ivBM-4PAP, no. 101098989), and Project ‘National Center for Gene Therapy and Drugs based on RNA Technology’ (CN00000041) financed by NextGeneration EU PNRR MUR—M4C2—Action 1.4—Call ‘Potenziamento strutture di ricerca e creazione di campioni nazionali di R&S’ (CUP J33C22001130001).