2.1.1. Describing crystal solvent channels with Voronoi diagrams

The problem of finding a path for the largest possible sphere is strongly related to the medial axis of the channel, which can be described as the subspace for which more than one closest point to the channel boundary exists (Blum, 1967). The medial axis includes all space where the center of maximal spheres can be placed, since a sphere touching only one boundary can still be enlarged in another direction.

In LifeSoaks, we use the closely related concept of a Voronoi diagram to describe the solvent-occupied space (Voronoi, 1908). For a collection of points, it can be shown that all corresponding Voronoi objects are also part of the medial axis (Fabbri et al., 2002). However, unlike the medial axis, a 3D Voronoi diagram can be calculated in expected in typical cases (Boissonnat & Attali, 2002; Amenta et al., 2007), even though the worst-case complexity is (Dewdney & Vranch, 1977). To compute the Voronoi diagram, we used an in-house implementation of a 3D incremental insertion algorithm to construct a Delaunay tetrahedralization (Green & Sibson, 1978; Guibas & Stolfi, 1985; Joe, 1991). A Delaunay tetrahedralization is dual to a Voronoi diagram. We can subsequently generate the Voronoi diagram in Θ(N) time from the tetrahedralization.

A Voronoi diagram in d dimensions computed on a set of input points P with assigns the subspace S_p to every point p for which the following relation holds:

Thus, S_p is the space for which p is the closest of the points in P. The resulting structure consists of convex polyhedra around every point p (see Fig. 1a for a 2D example). Their facets represent areas at which another point is equally distant, edges represent areas at which two additional points are equally distant and vertices are equally distant to four elements of P in total. The vertices and edges of the Voronoi diagram can be interpreted as a graph with an annotated radius of free space at every graph element. For vertices, this radius represents the distance to its four defining/closest elements in P. For edges, it is the minimal distance to its three defining elements and the radius is calculated from their planar circumcircle. The circumcircle also represents the intersection of the spheres of free space surrounding the vertices incident to the edge. The space occupied by all of these spheres represents our model of the solvent-occupied space.

Figure 1 (a) An example of a two-dimensional Voronoi diagram. The lines are parts of the perpendicular bisectors of the neighboring points and are good approximations for paths that maximize the minimal distance to points. (b) The same Voronoi diagram in a molecular interpretation. The points represent the atoms depicted in blue. The lines are paths with maximal distance to atom centers. The yellow line marks a potential channel with a certain minimum distance to the atoms. The edge radius is illustrated for edge e = (va, vb). re is the minimal distance of e to its closest atom center. For a molecular representation, re is reduced by ru, resulting in the minimal distance to an atom surface.

Constructing a Voronoi diagram from a molecular structure poses additional challenges in comparison to the general case. One of them is that we need to treat atoms as spheres for an adequate representation. However, since Voronoi diagrams in their standard form are generally defined on points, the available algorithms can only treat atoms as spheres of equal radius (see Fig. 1b). This limitation is sometimes addressed by approximating large atoms by several spheres (Yaffe et al., 2008; Chovancová et al., 2012) or even by sophisticated Euclidean Voronoi diagrams that treat arbitrary spheres explicitly, leading to curved Voronoi elements (Kim et al., 2005, 2015).

Since calculating Voronoi diagrams for a unit cell is computationally expensive, LifeSoaks uses a standard implementation with a unified sphere radius r_u, which will be subtracted from all vertex and edge radii. By default, r_u is set to the van der Waals radius of a carbon atom (1.7 Å; Bondi, 1964). While this only slightly overestimates the 1.5 Å radius of oxygen and nitrogen atoms, hydrogen atoms, with a van der Waals radius of 1.1 Å (Rowland & Taylor, 1996), are significantly smaller. Therefore, the corresponding spheres are moved r_u − 1.1 Å towards the connected heavy atom to decrease the overestimated space. In the common case that hydrogen atoms are not assigned in the PDB file, they will be added by our chemistry model with respective binding geometries beforehand (Urbaczek et al., 2013).

Another challenge is the common local coplanarity of substructures originating from symmetries or planar ring systems. Four atoms placed exactly on the same plane cannot form a tetrahedron, which may cause the algorithm to either fail or run for unfavorably long. This problem is addressed by perturbing the position of input points randomly within a 10⁻⁶ Å radius around their origin.

In an effort to reduce the computation time, the solvent-accessible surface (SAS) area is determined prior to calculation by the in-house implementation of a spherical-probe SAS algorithm (Richards, 1977; Reulecke et al., 2008). Atoms that are not solvent-exposed are discarded since they cannot be part of a solvent-channel boundary.

By default, small molecules and solvents such as water are not considered in the calculation and are treated as empty space. However, the user can include such molecules explicitly, for example to consider binders with high residence times. Like the protein itself, they are regarded as completely rigid since the model does not consider flexibility.

To represent entire crystal solvent channels by Voronoi elements, we need to construct a Voronoi diagram of one complete unit cell. Since a crystal can be constructed from one unit cell by applying translation operations, the same holds for its solvent channels. However, every element of the Voronoi diagram is defined by its closest atoms. Therefore, it is necessary to account for any atoms of neighboring unit cells that reach into the central unit cell or are close enough to the border to influence the outer Voronoi elements. This was achieved by adding additional symmetry copies at the unit-cell borders. Roughly half an additional unit cell in each direction and dimension is necessary to correctly compute the central unit cell, resulting in an input of approximately eight unit cells (see the supporting information for a detailed derivation of the space depending on the unit-cell geometry).

Since we are interested in a Voronoi diagram of one complete unit cell, the Voronoi diagram calculated from this input is subsequently cut at the borders of the central unit cell and artificial surface vertices are introduced at the cutting positions of edges. Furthermore, Voronoi elements are discarded if their radius is smaller than a user-defined cutoff since they are considered to be irrelevant for detecting solvent channels while inflating the resulting structure. Note that these artificial alterations result in a structure that is no longer an actual Voronoi diagram; therefore, the resulting graph will be referred to as Voronoi channel graph in the following.

2.1.2. Accounting for crystal periodicity in a Voronoi channel graph

In general, we are interested in calculating paths that traverse a macroscopic crystal consisting of numerous unit cells. Therefore, the paths that we calculate not only need to reach specific positions in the unit cell but also need to do so in a periodic manner. Starting from a crystal or unit-cell surface, the path, in theory, needs to reach infinitely many subsequent unit cells.

Such paths are trivial to detect in some cases, for example when the channel traverses a single unit cell and ends at its periodic starting point (Fig. 2a). However, as shown in Fig. 2(b), this condition is not always met. Even though the depicted channel traverses a complete unit cell, it is interrupted right after crossing the boundary to the next one and cannot reach its periodic starting point. These cases can be detected by projecting the intersections of opposing unit-cell surfaces onto each other and checking for overlap.

Figure 2 Four exemplary channel geometries that pose different requirements for the model. The channels are shown in yellow, while protein-occupied space is depicted in blue. (a) A trivial channel that linearly passes one unit cell before reaching its periodic starting point. (b) A channel that traverses the whole unit cell. Since its surfaces do not map onto each other, it does not traverse the macroscopic crystal. Therefore, it never reaches its starting point and cannot be detected by the model. (c) A channel that traverses the whole unit cell but crosses a second unit-cell boundary before reaching its starting position. These cases can be detected by considering the periodic boundary condition. (d) A special case that leads to the region marked in red not being correctly handled by the model.

In other cases, channels may traverse several unit cells before reaching their periodic starting point (Fig. 2c). Here, the Voronoi channel graph must accurately represent the periodic condition in the channel construction.

Therefore, we introduced artificial Voronoi edges that connect surface vertices directly in front of each other on opposite sides of the unit cell. On the macroscopic level, these edges represent connections to the neighboring unit cell. Subsequently, channels can be computed in full even if they cross the boundaries at some positions. A channel that exits the unit cell re-enters it on the opposite side when treated as passable.

A further special case in the model is shown in Fig. 2(d). It can occur when enclosed side areas of a channel are cut by the unit-cell boundary and these areas have no other connection to the main channel. In the algorithm, the unit-cell surface currently checked for channels has to be treated as non-passable, since otherwise the erroneous detection of local circles would be possible. Therefore, these special regions are isolated even though they are connected to the channel. In all cases visually inspected in the context of this study, this scenario had no significant impact.

2.1.3. Calculating bottleneck radii for channels

In general, the Voronoi channel graph may include a large number of paths from one crystal surface to the opposite surface. To predict the ability of molecules to traverse the crystal, we define the quality of each path by its narrowest subpart, i.e. its bottleneck radius. This decision leads to a particular difference from several previously mentioned Voronoi-based tools for detecting single protein channels, which often determine the optimal path using a shortest-path algorithm. Their distance measure often includes an inverse contribution of the local radius. However, the resulting optimal path also depends on its length, which we wanted to avoid.

Instead of computing the shortest paths, we define channels with respect to an arbitrary radius r′ as a subset of Voronoi vertices that can be reached by a sphere with radius r′ on an infinite periodic path through the crystal.

Let I_P be the set of all infinite crystal paths that periodically traverse subsequent unit cells, thereby describing a path through the macroscopic crystal. For a given vertex v the following relation holds:

Thus v is reachable by a sphere of radius r′ if and only if there exists an infinite path in the Voronoi channel graph that contains v while the radius r(e) of its thinnest edge is at least r′.

The LifeSoaks algorithm aims to annotate each vertex with its local bottleneck radius r_b(v), for which the following relation holds:

Thus r_b(v) is the largest bottleneck radius of all paths that include v. This radius may be much lower than the local vertex radius r(v), for example if v is part of a large enclosed cavity that can only be reached by a small channel. However, r_b(v) can never exceed r(v).

r_b(v) can be stored locally at each vertex. A channel set is derived by including every vertex with r_b(v) ≥ r′. The principle is visualized in Fig. 3.

Figure 3 A two-dimensional schematic of a crystal unit cell with a simplified Voronoi channel graph. Blue objects are proteins with a ligand in a binding site. The double arrows mark bottleneck positions. Bottleneck a with radius ra defines the space visualized in dark yellow that a sphere of radius ra or smaller can reach. Bottleneck b has radius rb. It defines the space marked in light yellow. Since rb < ra, a sphere of radius rb can also be placed in all of the dark yellow space. The white space cannot be reached by either of the two spheres but may potentially be reached by smaller spheres. A typical unit-cell Voronoi diagram consists of thousands of edges and vertices. The exact number can optionally be reported in the output.

This concept is realized by implementing an algorithm that is based on a union–find data structure that saves disjoint sets of Voronoi vertices in a new directed tree graph. Finding the set of a vertex as well as combining two sets into one can be efficiently performed in , where m is the number of set elements and α denotes the inverse of the Ackermann function, which can be considered to be constant for practical inputs (Ackermann, 1928; Galler & Fisher, 1964; Tarjan & van Leeuwen, 1984). The pseudocode is a simplification of the main routine, neglecting several implementation details such as optimizations. A schematic visualization of the main steps of the algorithm can be found in Fig. 4.

Figure 4 A two-dimensional schematic of the LifeSoaks algorithm. The unit cell includes two internal vertices, four surface vertices, five internal edges and two artificial surface-connecting edges marked *. We note that and due to internal point symmetry. Further, r(e1) > r(e2) > r(e3). (a) The state before any edge is processed. In this iteration the bottleneck radii are determined from left to right, so vertices on these sides are marked as surface vertices (magenta coloring) and e*1 is excluded. (b) After two iterations e1 and have been processed, since they are the largest edges. Because they connect one new vertex to a surface vertex, their incident vertex is now also accessible by one surface but not by both. (c) After three additional iterations e2, and e*2 have been processed. Since e*2 connects the surface vertices at the top and bottom, all vertices are now part of the same set. Since this set now contains vertices of the opposing surfaces, it represents a unit-cell-traversing path and all vertices will be annotated with the bottleneck radius r(e2). (d) Finally, e3 is processed. However, because the incident vertices are already part of the same set, no set or radius is changed.

The algorithm is performed three times, once for every unit-cell dimension/axis (line 3 in the pseudocode shown below). In each iteration, one pair of opposing unit-cell surfaces is checked for paths through the entire unit cell. Therefore, artificial surface-connecting edges are not included in the calculation for this surface pair (line 7). The other two surface plane pairs are treated as passable to detect channels, as depicted in Fig. 2(c).

To begin with, every Voronoi vertex is a single set. Surface Voronoi vertices of the two considered planes are labeled accordingly. Afterwards, all Voronoi edges are processed in descending order with respect to their radii (lines 2 and 5). Whenever an edge connects two previously disjoint union–find sets, these sets are combined (line 9) and whether this set contains a path through the unit cell is checked. This holds if one of two conditions is met. If one set already includes a complete path, this path is now accessible to the union of both sets. However, if neither set includes a previously discovered path, it is checked whether they reach opposing sides of the considered surface pair and whether these surface vertices can be projected onto each other (line 16). In both cases the union represents an infinite path through the crystal and the radius of the current edge is annotated to every included Voronoi vertex unless it had previously been assigned a larger radius (lines 17–20). The radius-descendent edge processing ensures that a connected path is always found exactly when its narrowest edge is included, automatically delivering its bottleneck radius for that dimension.

Ultimately, every vertex is annotated with a radius for each unit-cell dimension, the largest of which represents its bottleneck radius. For a specific vertex v, this radius r_b(v) corresponds to the radius of the largest spherical object that can reach this vertex on an infinite periodic path through the crystal. The largest local bottleneck radius among all vertices corresponds to the bottleneck radius of the main channel of the crystal. Note that separate bottleneck information for all three dimensions can also be of interest since the blocking of diffusion directions can be associated with significantly increased soaking times (Geremia et al., 2006).

2.1.4. Assigning bottleneck radii to binding pockets

To describe the solvent-channel path with respect to a binding pocket, we introduced a routine that assigns a bottleneck radius to user-defined binding sites.

In our code base, binding sites are either derived from a unique ligand identifier by defining every residue within a 6.5 Å radius as part of the pocket or directly provided by additional PDB files containing the binding-site residues. The convex hull for all atoms of pocket residues is calculated and defined as the pocket space.

If the pocket space contains Voronoi elements, a local bottleneck radius could be extracted from this information alone. However, the complex binding event is not well described by our simplified representation. Instead, we want to predict whether the ligand can reach the space in front of a binding pocket. Therefore, the Voronoi channel graph is searched for the largest bottleneck radius r_b,poc within a user-defined distance of the convex hull of the pocket. After determining r_b,poc, a breadth-first search is performed on the Voronoi channel graph to locate the corresponding bottleneck position, i.e. the first edge e = (u, v) for which r(e) = r_b,poc = r_b(u). Note that this method may fail if the binding site is not pocket-like since it might contain no Voronoi elements. Observed failures can usually be attributed to flat convex hulls of binding sites on the protein surface.

2.1.5. Data representation and visualization with LifeSoaks

LifeSoaks assigns a bottleneck radius to every Voronoi vertex and to the crystal structure as a whole. To allow dynamic visualization of channels considering all local bottleneck radii, LifeSoaks writes grid files in the CCP4 format (Winn et al., 2011) by converting the annotated Voronoi channel graph into a grid. Each grid point is annotated by the local bottleneck radius of the largest vertex sphere containing it. Various molecular-visualization programs such as UCSF Chimera (Pettersen et al., 2004), PyMOL (Schrödinger) and Coot (Emsley et al., 2010) support this format and display the grid data with respect to user-defined cutoff values. Note that for accurate channel representation the true values have to be visualized and not the standard deviation σ of the grid, which is commonly used for electron-density maps. For a specific radius, the displayed grid volume represents the space in which a sphere of that radius could roll around freely on an infinite path through the crystal. The radius represents r′ and the visualized grid resembles the corresponding as described in formula (2).

With respect to a position of interest, such as a binding pocket, the user may gradually decrease r′ until the channel reaches that position.

In addition to the standalone command-line tool, the functionality is freely available on the ProteinsPlus web server at https://proteins.plus (Fährrolfes et al., 2017; Schöning-Stierand et al., 2022).

Starting with LifeSoaks version 1.1, the command-line tool will also include an option for two additional types of grid visualization. One contains the true local radius at each position, so a sphere with the current cutoff radius may be placed anywhere inside the grid but can potentially not move around freely. The other contains the local distance to the closest atoms. The displayed grid represents the space in which the center of the potentially blocked sphere can be placed.

2.1.6. Matthews coefficient and solvent-content calculation

To compare the bottleneck radii with the crystal solvent content and the Matthews coefficient, we did not use the values deposited in the PDB since we discovered several inconsistencies in the deposited values. We suspect that slight differences in the various software tools used for crystallo­graphic structure determination are the cause of this. We recalculated these values from the protein mass M_a of our internal chemistry model (Urbaczek et al., 2011, 2012) to obtain consistent results. The solvent content is calculated using an average protein density δ of 0.81 Da Å⁻³ as described by Matthews (1968). It is then expected to be the volume not occupied by protein p_solv = 1 − p_prot = 1 − [M_prot/(δ × V_a)]. V_a represents the volume per asymmetric unit. Note that, for the sake of simplicity, the calculation is not altered for nucleic acids, which are generally assumed to have a lower density (Kantardjieff & Rupp, 2003).

2.1.7. Implementation and computational requirements

LifeSoaks was written in C++ as part of the NAOMI cheminformatics library (Urbaczek et al., 2011, 2012). The computation of the Delaunay tetrahedralization, which is not parallelized, takes the major part of the computation time and storage. The conversion into a Voronoi channel graph and the union–find operations of the channel-computation routine can be neglected in comparison. The tool typically requires a few minutes for one crystal structure on a standard desktop machine. Longer computation times occasionally occur for large structures with complex crystal symmetries or due to unfavorable input distributions. Cases with Ω(NlogN) tetrahedra also lead to high memory requirements. LifeSoaks has an included heuristic memory limitation which terminates the calculation at a user-defined memory usage to prevent memory overflows. This is set to 6 GB by default but needs to be increased for large structures.

2.2.1. Data-set generation

We designed a data set of protein–ligand complex structures that were obtained by soaking a ligand-free crystal in a solution containing the ligand. To this end, we screened the REMARK 280 lines of all X-ray crystallographic structures in the PDB. These remark lines contain experimental details of the crystallization conditions. Structures containing the string `soak' (case-insensitive) in this remark were analyzed. Based on the description of the crystallization procedure in the PDB file, we built a data set including the PDB codes and the HET codes of the soaked ligands. We extracted the experimental details from the associated publication whenever the description was ambiguous. Finally, we split the data set into two subsets: (i) soaked small molecules, ions and atoms (`Ligand-Soaked Complexes') and (ii) soaked oligomeric or polymeric compounds, i.e. oligosaccharides, peptides, proteins or nucleic acids (`Oligomer-Soaked Complexes'). The second subset was not used in the subsequent analyses as the reliability of conformer predictions for oligomeric or even polymeric molecules with many rotatable bonds has limitations (Friedrich et al., 2017), leading to poor interpretability of the results. The first subset used in our analyses comprised 1713 structures and 4718 ligands that have been reported to have been soaked into preformed crystals. Whenever these ligands were found in the structure, we considered all corresponding ligands and their binding sites. For ligand detection, we used an in-house tool developed with our NAOMI library (Urbaczek et al., 2011, 2012). All entities identified as ligand molecules and residues listed in the HET section of the PDB files are checked for the HET code of the soaked molecule. The ligand instances with this HET code were considered in soaking predictions. The full data set, including the PDB codes and ligand identifiers, is readily available via the GitHub repository for the project (see the Supporting information).

2.2.2. Ligand radii assignment for ligand-soaked complexes

Since our geometrical model describes the objects traversing a channel by a single radius, we need a procedure to determine a meaningful radius for a given molecule conformation. We first predict the most probable protonation state of each ligand from its SMILES representation with UNICON (Sommer et al., 2016; Sommer, 2016). The protonation and tautomer-generation mode was set to single, i.e. only the top-scoring states are used. Subsequently, the conformer-generation method Conformator (Friedrich et al., 2019; Friedrich, 2019) was applied to generate at most 1000 conformations per ligand with default parameters. Both tools are freely available for noncommercial and academic use and can be found using the links in the corresponding references. For seven ligands (PDB ligand IDs AF3, BEF, DQY, PA0, SEY, VM4 and WO4), preparation of the ligand with UNICON failed and only their ideal conformation as downloaded from the PDB was used. However, for one ligand (ligand ID 6VQ) the workflow was successfully performed with the respective input in SDF format with ideal coordinates from the PDB but not with SMILES.

Based on the number of generated conformations, we split the set into molecules with a single conformer and molecules with multiple conformers. The minimum projection radii of all conformers were subsequently calculated using the ChemAxon tool cxcalc (https://www.chemaxon.com). Based on the van der Waals radii of the molecule atoms, the method searches for the enclosing circular planes onto which all atoms can be projected. It reports the radius of that with the smallest radius (the minimum projection radius). For single atoms (for example xenon) and ions that could not be annotated with their respective projection radii, we referred to two publications reporting experimentally determined and calculated radii (Shannon, 1976; Rahm et al., 2016). The assigned projection radii can be found in the GitHub repository for the project.

2.2.3. Software parameters

LifeSoaks runs on the hand-curated benchmark data set were performed with the default settings apart from the following changes: the minimum channel radius was adjusted to 0.5 Å if the minimum projection radius of a ligand was below the default of 1.7 Å and the memory limit was set to 120 GB. None of the calculations for the data set exceeded this limit. Protein structures in PDB format were used to reproduce the behavior of the tool on the ProteinsPlus web server. PDB structures with multiple models lead to software failures and should be preprocessed by the user to include only one model. MAP_CHANNELS runs were performed with the default settings apart from the following changes: we enforced the use of the specified grid size of 2.0 Å regardless of the maximum time. In addition, we set the channel cutoff, i.e. the minimum grid–atom surface distance to be considered a channel, to 0.5 Å.

2.2.4. EDIA calculations

Since we work with a rigid protein model but expect that highly flexible residues do not represent a strict obstacle for soaked molecules, we performed a second calculation in which we excluded parts of the structure that are poorly supported by electron density. To this end, the electron-density score for individual atoms that quantifies the electron-density fit of an atom was calculated with EDIA (Meyder et al., 2017). If available, the corresponding 2F_o − F_c electron-density maps were downloaded from the PDBe (Armstrong et al., 2020). For structures with available electron-density maps, EDIA was applied with default parameters. Residues with at least four atoms with an EDIA score below 0.8 were excluded from the solvent-channel calculations with LifeSoaks in this additional run.

3.1.1. Radius distribution and correlation with the Matthews coefficient

To estimate the benefit of predicting a solvent bottleneck radius for a user, we compared this value for every considered structure with the corresponding solvent content and Matthews coefficient (Fig. 6). Despite the inverse relation between the Matthews coefficient and the solvent content, both comparisons are instructive. The Matthews coefficient, on the one hand, has the advantage that structures with a large solvent content are better distinguished from each other since the volume per mass has no upper limit. The solvent content, on the other hand, is more intuitively interpretable.

Figure 6 The correlation between the overall bottleneck radius and the solvent content (a) and the Matthews coefficient (b) in the form of scatter plots with univariate histograms. Outlier values with a Matthews coefficient greater than 8 Å3 Da−1 or a bottleneck radius greater than 80 Å are not displayed for improved visualization.

Both parameters correlate with the bottleneck radius, illustrating that a structure with high solvent content is likelier to have large channels. In particular, large channels with bottleneck radii of more than 20 Å can almost exclusively be found in structures with greater than 50% solvent content. However, the relation is far from strict, as supported by the linear Pearson correlation coefficient of 0.61 between the bottleneck radius and the Matthews coefficient. In general, many examples of narrow bottlenecks paired with comparably large Matthews coefficients can be found. At the same time, other structures show tubular main channels in a heavily protein-occupied unit cell, leading to a low Matthews coefficient but a large bottleneck radius.

3.1.2. The example of phosphoribulokinase from Rhodobacter sphaeroides

An example of a very small bottleneck radius paired with a very large solvent content was observed in the structure of phosphoribulokinase from R. sphaeroides crystallized in the cubic space group P432 (Harrison et al., 1998; PDB entry 1a7j). Considering the solvent content of 57.2% and the Matthews coefficient of 2.88 ų Da⁻¹, one would initially not expect an obstacle for soaking. However, the bottleneck radius of the main channel is only 2.24 Å. Since this is an extremely small radius for a crystal with such a high solvent content, we inspected the structure visually. As shown in Fig. 7, the 24 symmetry-related copies of the protein form a spherical enclosure of the solvent-occupied space in the middle of the unit cell. This space, however, is only reachable by a very narrow path from each unit-cell surface, which exhibits a bottleneck defined by four symmetric α-helices from four protein copies. The closest residue in the α-helices is Ser248.

Figure 7 Cubic crystal of phosphoribulokinase from R. sphaeroides (PDB entry 1a7j) with its channel volume visualized with a cutoff radius of 2.24 Å. (a) The whole unit cell with the main channel. The 24 symmetric copies enclose a sphere of solvent that is only reachable by a small channel with a bottleneck radius of 2.24 Å. (b) The corresponding bottleneck, in which two of the four Ser248 residues are depicted. Images were created with UCSF Chimera (Pettersen et al., 2004).

3.2.1. Predicting the soakability of binding sites with single-conformer ligands and ions

For 98.5% of the binding sites, we could correctly predict their accessibility to small molecules when only considering the minimum projection radii, with only 30 sites not predicted to be soakable. In eight cases, the main channel bottleneck radius was not predicted to be large enough to allow the ligands to traverse the crystal (Fig. 8a). These eight human oxidized purine nucleoside tri­phosphate hydrolase structures were determined in complex with the fragment 1H-imidazo[4,5-b]pyridin-2-amine (PDB ligand ID BU8) with a minimum projection radius of 3.19 Å. According to the experimental data section of the PDB files, the fragment was soaked into the crystal. However, a publication is lacking for these examples, so we could not further assess whether this agrees with the actual experimental conditions. A closer examination of the solvent channels in the structure of this enzyme (PDB entry 6glg) shows that the bottleneck is formed by poorly resolved loop regions that might be sufficiently flexible to allow an increase in the bottleneck radius. These regions spanning residues 24–28 and 141–143 are characterized by high B factors, indicating a flexible region. Its residue movements might enlarge the solvent channel. Indeed, the bottleneck radius of the main channel of the crystal and that in front of the binding site increase from 2.98 Å to 4.12 and 3.95 Å, respectively, if residues with at least four atoms with an EDIA (Meyder et al., 2017) below 0.8 are not considered in the LifeSoaks calculation (Fig. 8c).

Figure 8 Crystal and binding-site bottleneck radii and their relation to the minimum projection radii of molecules with a single calculated conformer. (a) The crystal bottleneck radii are plotted against the minimum projection radius of the individual soaked molecules. (b) The bottleneck radii observed in front of the small-molecule-binding site are plotted against the minimum projection radius of the molecules. The line of identity is provided to highlight false-negative results. The points are colored according to the kernel-density estimate using Gaussian kernels visualizing highly occupied regions in the data set. (c) Loop flexibility resulting in false-negative predictions with LifeSoaks as shown for the structure of human oxidized purine nucleoside triphosphate hydrolase (PDB entry 6glg). The atom and cartoon representations are colored according to their B factor and the black sphere indicates the bottleneck position. The channels predicted by LifeSoaks are depicted in yellow with a cutoff radius of 2 Å. (d) Insufficient electron-density support for bottleneck-constituting residues leading to false-negative LifeSoaks predictions, as illustrated using the structure of mevalonate diphosphate decarboxylase (PDB entry 6e2v). The electron density (2Fo − Fc) is shown as a grid at a sigma level of 1. Residues within a 5 Å environment of the bottleneck position (black sphere) are shown as sticks and their atoms are colored according to their B factor. (c) and (d) were generated with PyMOL (Schrödinger).

For the remaining 1955 structures, we analyzed the accessibility of the small-molecule-binding sites from channels in the crystal (Fig. 8b). We visually inspected the structures of the 22 false-negative predictions. In general, the differences between the minimum projection radius of the ligands and the bottleneck radius of the corresponding site are very small. They can most likely be attributed to small channels leading to the binding site, which is on the edge of our simplified model of a ligand diffusing through the crystal, as passage also depends on the local radii of the ligand. The difference is greater than 0.5 Å in only eight cases, which will be discussed below.

In the structure of leukotriene A4 hydrolase with PDB code 3fuh (Davies et al., 2009), the ligand with the PDB three-letter code 5H1 (1H-indol-5-ol, minimum projection radius of 3.69 Å) is nearly fully enclosed by the protein chain, with only a narrow channel leading to its binding site. A loop region closes over the ligand and the loop residues might be flexible, enabling the ligand to reach the binding site in a more open conformation. Moreover, the enzyme has high intrinsic flexibility and is characterized by open and closed states (Stsiapanava et al., 2017). These conformational differences cannot be captured only considering the rigid crystal structure.

In another case, one sulfite ion (, PDB three-letter code SO3 in chain A with residue sequence number 537) of the cytochrome c nitrite reductase structure with PDB code 3lg1 (Trofimov et al., 2012) is flanked by the helices of three asymmetric units and the bottleneck is formed by Val454 in chain A. Intriguingly, the helix ends constituting the bottleneck are characterized by high average residue B factors. A similar situation is observed for the structure of the same enzyme (PDB entry 2zo5; Polyakov et al., 2009) and the corresponding binding site of a nitrite ion (PDB ligand code AZI with residue sequence number 530 in both chains A and B). The termini of the helices closing the binding site are characterized by high temperature factors. This flexibility might enable entrance of the ion into the binding site through a channel.

Similar observations were made for several structures of cytochrome c peroxidase [PDB entries 4jm8, 4jm6, 4jma, 4jmw (Rocklin et al., 2013), 1aed, 1aeo, 1aeg, 1aem, 1aee, 1aen, 1aef, 1aeb (Musah et al., 2002), 1aeu, 1aet and 1aes (Fitzgerald et al., 1996)] with PDB ligand IDs 26D, LG3, DTI, 3FA, 2AP, IPH, 4AP, MPI, ANL, 25T, 2MZ, 3AP, 1MZ, 3MT and IMD, respectively. A loop region spanning residues 190–195 encloses the binding site of the respective ligands. This loop was reported to be highly flexible in the apo structure of the protein (Fitzgerald et al., 1996). For the PDB entry 4jm8 (Rocklin et al., 2013), we observe a bottleneck with a radius of 2.92 Å, which is smaller than the minimum projection radius of the ligand (3.87 Å). In addition, this bottleneck is not situated at the entrance to the corresponding site but in a channel blocked by heme. After omitting residues with low EDIAm values, which indicate poorly defined electron density and high flexibility, the bottleneck radius increases to 9.89 Å. Thus, the pocket is reachable by the solvent channels in the crystal. In the ligand-occupied state, this loop closes and the ligand is located in the center of a very narrow channel.

For two caesium ions (ligands with residue name CS and residue sequence numbers 502 and 501 in chains B and C, respectively) in the structure with PDB code 6hfb (Zakrzewska et al., 2019), the bottleneck radii in front of the binding site were calculated to be smaller than their radii. However, the radii inside and in front of the pocket are identical, explaining this finding: these two ions lie on the protein surface and their sites are defined by only two residues, rendering the binding-site determination by LifeSoaks unreliable. However, we do not consider this to be a limitation of the method because users are expected to be mainly interested in druggable and sufficiently buried binding sites.

Examples of the latter four typical challenges for predicting the crystal soakability and binding-site accessibility are depicted in Supplementary Fig. S4.

The last case is the binding site of a beryllium trifluoride ion (minimum projection radius of 2.66 Å) in a structure of mevalonate diphosphate decarboxylase (PDB entry 6e2v; Chen et al., 2020). The side chain of Lys67 with very high atomic B factors forms the bottleneck in front of the pocket (2.65 Å). Its side-chain conformation lacks proper electron-density support and might be present in various conformations in the crystal (Fig. 8d).

In summary, LifeSoaks correctly predicted crystal soakability and binding-site accessibility for most of the true-positive examples. However, if the proteins under investigation contain highly flexible regions, the predictions should be carefully evaluated. In particular, flexible protein regions that might block the binding site in the ligand-occupied state but are more flexible in the ligand-free structure must be considered when assessing the performance of the tool. All cases in which our tool fails to correctly predict the soakability are connected to intrinsic features of the corresponding protein structure models, such as flexibility and poor electron-density support for atoms. The latter might be attributed to high flexibility or an overall poor crystal structure resolution.

3.2.2. Predicting the soakability of binding sites with multiple-conformer ligands

A slightly different scenario was observed for larger molecules with multiple conformations (2445 ligands; Figs. 9a and 9b). Bottleneck radii inside and outside the pocket could not be calculated for six sites. In 46 cases, the binding-site-defining point was predicted to be inaccessible (completely enclosed by the protein). Therefore, a bottleneck outside the pocket was not calculated. For the remaining 2393 structures in this subset, the calculated bottleneck radii (crystal and binding-site bottleneck radii) are often smaller than the minimum projection radius of the conformational ensemble of the corresponding molecules compared with the set of rigid small molecules (231 false-negative predictions). For 148 of these false-negative predictions, the bottleneck radius of the crystal was predicted to be smaller than the minimum projection radius of the conformational ensemble of the soaked ligand. However, most data-set structures were correctly predicted as soakable (90.1%).

Figure 9 Crystal and binding-site bottleneck radii and their relation to the minimum projection radii of the conformational ensembles of the soaked small molecules. (a) The crystal bottleneck radii are plotted against the minimum of the minimum projection radius of the conformational ensemble of the individual soaked molecules. (b) The bottleneck radii observed in front of the small-molecule-binding site are plotted against the minimum of the minimum projection radius of the conformational ensemble of the molecule. The line of identity is provided to highlight false-negative results. The points are colored according to the kernel-density estimate using Gaussian kernels visualizing highly occupied regions in the data set. (c) The modeled structures of flexible loops that lack proper electron-density support, leading to false-negative predictions with LifeSoaks as shown for the structure of the human protein kinase CDK2 (PDB entry 4ez7). The atom and cartoon representations are colored according to their EDIAm values. The electron density (2Fo − Fc) is shown in grid representation at a sigma level of 2. (d) Huge domain movements because of induced-fit phenomena resulting in false-negative LifeSoaks predictions, as illustrated by the structure of prolyl endopeptidase from Aeromonas caviae (PDB entry 3muo, blue). The original structure of the proteins in the crystals used for soaking highlights the accessibility of the binding site in the ligand-free crystal of the structure (PDB entry 3ium, cyan). The binding site is highlighted by the corresponding ligands in stick representation (ZPR_A_701 and ZPR_B_702, magenta). (c) and (d) were generated with PyMOL (Schrödinger).

The higher number of structures whose soakability was not correctly predicted may be attributed to the molecules being larger than those in the first data set. This finding highlights the boundaries of the basic assumption that we can predict soakability by approximating the minimum dimensions of a molecule by a 2D-projected circle of minimum radius and annotating channels by spheres. Additionally, inaccuracies in the conformational sampling, especially in the case of large molecules with many rotatable single bonds, might occur. Finally, some of the proteins whose flexibility was discussed earlier can also be found in this data set, only with different ligands (for example human oxidized purine nucleoside triphosphate hydrolase structures, PDB entries 6gle, 6gli, 6gln, 6glr, 6gls and 6glv). Others are new challenges for soaking prediction that can only be faced by careful visual inspection of the protein structures under investigation. To this end, we sorted the results according to the difference between the bottleneck and the minimum projection radii in descending order. An analysis was performed starting with the prediction with the highest difference between the minimal projection radius and the predicted bottleneck radius. We will discuss some examples with especially large differences.

The example with the highest difference between the bottleneck and the minimum projection radius is the protein kinase structure with PDB code 4ez7 (Martin et al., 2012). Both the crystal bottleneck radius and the bottleneck radius of the channel leading to the binding site for 8-anilinonaphthalene-1-sulfonic acid (residue 2AN with residue sequence numbers 302 and 303 in chain A) and staurosporine (residue STU with residue sequence number 301 in chain A) are smaller (2.16 Å) than the minimum projection radii of the soaked molecules (4.62 and 5.63 Å, respectively). However, large regions of the modeled structures are not properly supported by electron density due to low resolution and high flexibility (Fig. 9c). The arrangement of the asymmetric units in the unit cell might also allow loop flexibility. The ligand is flanked by loop regions with high temperature factors and no reliable electron density can be found for the modeled activation loop of the kinase. Upon removing un­reliably modeled residues according to their EDIA value (at least four atoms with EDIA < 0.8), the bottleneck radius of the structure increases to 6.55 Å, which can be attributed to the omitted loop regions opening a broader channel through the crystal (Fig. 9c). These results indicate that low electron-density support causes reduced prediction reliability. The same holds for the protein kinase structures with PDB codes 3pxz, 3py1 and 3pxq (Betzi et al., 2011).

Another interesting example is the structure of (S)-1-phenylethanol dehydrogenase (PDB entry 2ewm; Höffken et al., 2006). It has a very small crystal bottleneck radius of 2.63 Å. Intriguingly, a structural region spanning residues 188–206 of chain B is characterized by very high B factors. Unfortunately, an electron-density map is not available from the PDBe (Armstrong et al., 2020) to check the validity of the modeled atom coordinates. However, a ligand-free structure of this enzyme (PDB entry 2ew8) crystallized in a different space group and is characterized by different unit-cell parameters. In this structure, residues 188–206 are not resolved, leading to a much higher calculated crystal bottleneck radius of 9.5 Å. According to the authors of the structure (Höffken et al., 2006), the crystal used to obtain the NAD-bound structure was grown under the same conditions. It broke apart upon soaking and one piece was immediately cooled for structure determination. Therefore, it cannot be excluded that the ligand traversed the crystal in the apo form with the former space group, giving rise to changes in the crystal after soaking. Upon omitting the discussed flexible regions from the LifeSoaks calculations, the calculated bottleneck radius increases to 4.0 Å, which is still smaller than the minimum projection radius of the soaked NAD molecule (5.04 Å). Omitting another region that was not resolved in the ligand-free structure and was characterized by high B factors in the soaked structure (residues 188–192 of chain A) leads to a bottleneck radius of 4.4 Å that is at least similar to the predicted minimum projection radius of NAD.

The structure of prolyl endopeptidase (PDB entry 3muo) exhibits a large crystal channel of 22.8 Å. Nevertheless, its binding sites are inaccessible due to high enclosure by the protein. Thus, there must be another reason for its soakability than local protein flexibility. In the crystallized opened state of the enzyme (for example PDB entry 3ium; Li et al., 2010) its two domains are apart from each other, leading to highly accessible binding sites (Fig. 9d). These open-state crystals were used to soak ligands into the structure (Li et al., 2010), although the soaked structure is characterized by different crystal properties. Indeed, the authors of the structure explained this observation. They argue that due to the high solvent content of the crystal (62%), huge conformational changes of the enzyme are well tolerated. They do not interfere with crystal integrity due to the loose packing. These differences explain why LifeSoaks could not correctly predict the soakability of the crystal in the ligand-bound form. Therefore, it is advisable to exploit structures of the native crystal to be soaked for predictions with LifeSoaks. The same holds for the structure of this enzyme for PDB entry 3ivm (Li et al., 2010; 57% solvent content).

Similarly, we find a large solvent channel that traverses the crystal in the structure of aminopeptidase N with PDB code 2dqm (Ito et al., 2006). However, the binding site of the soaked small molecule bestatin is connected to this channel by a narrow channel with a bottleneck radius of only 1.8 Å, which is much smaller than the minimum projection radius of the soaked molecule (3.95 Å, residue name BES). The bottleneck is formed by Lys852 and Arg845 from two neighboring asymmetric units. These residues are characterized by higher B factors than the average for the structure and might be sufficiently flexible to allow ligands to traverse. Considering this flexibility and the high solvent content of 66%, it might be possible that conformational changes accompany the binding of the ligand. Different conformational states of aminopeptidase N are also discussed in the literature (Addlagatta et al., 2006).

The structure of lysozyme in complex with the soaked dye bromophenol blue (PDB entries 6syc, 6syd and 6sye, tetragonal space group; Plaza-Garrido et al., 2020) also represents an interesting explanation for false-negative predictions regarding the ability of small molecules to traverse a crystal. The crystal bottleneck radii are considerably smaller than the minimum projection radius of the dye (5.4, 4.8 and 3.8 Å, respectively, versus 5.8 Å). The bottleneck radii in front of the binding site are even smaller, contradicting the explanation in the corresponding publication (Plaza-Garrido et al., 2020). Its authors used MAP_CHANNELS and the radius of gyration of the dye (4.6 Å), which underestimates the molecular dimensions, to analyze the accessibility of two sites in the crystal. However, our calculations cannot confirm that the soakability can be explained without considering the flexibility of the protein. Intriguingly, a larger bottleneck radius of 5.4 Å was found for the structure in an orthorhombic space group (PDB entry 6syc). For the corresponding ligand-free orthorhombic (PDB entry 6f1o) and tetragonal (PDB entry 6f1p) lysozyme crystals (Plaza-Garrido et al., 2018), similar bottleneck radii were predicted. As the bottleneck is lined by less flexible proline and serine residues (Pro79 and Ser86), there might be another explanation for the soakability of this dye molecule. The nonspherical 3D conformation of the soaked molecule stands out. The assumption that a single minimum projection radius describes a molecular conformation is too simple to explain the soakability of the crystal. An ellipsoid describes the form of the molecule much better and might lead to better predictions. Also, different conformations of the molecule might play a crucial role for a ligand traversing a crystal. Although we analyzed the minimum projection radii of the conformers to obtain the smallest one, we cannot exclude that a more comprehensive conformer sampling might lead to conformers with even smaller minimal projection radii. Also, local minimal projection radii for various conformations might play a role in crystal traversal of small non-globular-shaped molecules.

In summary, the reasons why LifeSoaks might fail to predict that crystals can be soaked with ligands are high local flexibility, domain movements, missing electron-density support for modeled structures, a missing ability to cope with PDB structures containing multiple models, an insufficient sampling of the small-molecule conformations leading to inaccurately predicted projection radii, an inadequate representation of the solvent space and ligands by spheres, and the inability to detect completely enclosed sites, for example due to induced-fit phenomena. With regard to the former issues, users might modify the corresponding structures accordingly by excluding highly flexible residues (based on the B-factor distribution or EDIA calculations) to exclude unreliably modeled regions based on the electron density. Furthermore, the use of structures of the protein in the ligand-free crystals that are to be soaked is advisable. Also, users should carefully inspect the predicted channels and visualize the bottleneck positions to assess the true soakability of their crystals manually. The remaining issues cannot be fully solved without substantial further developments. However, handling protein flexibility is one of the most important issues to ensure reliable predictions in most analyzed examples.

3.2.3. Comparing LifeSoaks and MAP_CHANNELS

To the best of our knowledge, MAP_CHANNELS (Juers & Ruffin, 2014) is the only other available software for automated bottleneck calculations and soakability predictions. Although analyses of binding-site accessibilities are not feasible with MAP_CHANNELS, we decided to compare both tools with regard to the calculation of crystal bottleneck radii for the PDB structures in the data set. Fig. 10 shows the crystal bottleneck radii of the set of small-molecule-soaked crystal structures as predicted by LifeSoaks and MAP_CHANNELS. Major differences can only be observed in the 2D and 3D radii of the structures (Supplementary Fig. S5). The 2D and 3D radii describe the largest radius that a sphere may have to move infinitely along the 2D and 3D directions, respectively, in the crystal. They are smaller than or equal to the overall bottleneck radius, which can also be called the 1D radius. With regard to the 1D radii (the bottleneck for traversing the crystal in one direction), which are the most informative radii for a user investigating the general soakability of a crystal, both tools predict very similar bottleneck radii. In nine cases, the bottleneck radii are predicted to be more than 2 Å smaller by LifeSoaks. However, these cases represent crystals with very large bottleneck radii of 13 Å and above. In only one case, MAP_CHANNELS predicts a considerably smaller bottleneck radius (>2 Å difference). This is the structure of decay-accelerating factor (PDB entry 1ok9; Lukacik et al., 2004). The outstanding feature of this structure is the expansion of the asymmetric unit through multiple unit cells, which might lead to inaccuracies in the channel analyses with LifeSoaks, leading to a difference of 2.1 Å.

Figure 10 Comparison of the 1D bottleneck radii calculated by LifeSoaks and MAP_CHANNELS (1713 structures). The plot covers most data points up to 20 Å. The blue dashed lines indicate a prediction difference of at most 2 Å.

In general, MAP_CHANNELS analyses lead to lower bottleneck radii compared with LifeSoaks, within the limitations of the chosen grid spacing for MAP_CHANNELS calculations. The resulting maximal radius underestimation is half of the maximal distance that any point in space can have to its closest grid point. This extreme case would be equivalent to a point exactly in the center of an eight grid-point cube, resulting in an error of ( Å for our analysis). Another intriguing difference between the methods is the treatment of protein modifications, as exemplified by the predicted bottleneck radii for the structure of chlorperoxidase from Leptoxy­phium fumago (PDB entry 2ciz). The protein has several glycosylation sites, leading to a bottleneck radius of 10.5 Å as predicted by LifeSoaks. In contrast, MAP_CHANNELS predicted a radius of 12.1 Å, neglecting the glycosylation of the protein. After re-performing the calculation enabling the option to include nonwater heteroatoms in the distance-map calculation, the method reports a similar radius of 10.6 Å. However, a user would expect that the calculation would consider all covalent protein modifications as part of the protein. Additionally, explicitly using all heteroatoms apart from those of water molecules might lead to problems due to considering additional buffer components that can freely diffuse in the crystal in the calculation.

For the complete data set, the MAP_CHANNELS crystal bottleneck radii of 128 structures suggest nonsoakability of the crystals, compared with LifeSoaks, which predicted only 103 structures to be nonsoakable based on the crystal bottleneck radii. Notably, the accuracy of MAP_CHANNELS can be increased with a change in grid spacing at the expense of a longer computing time.

Regarding run time, the most significant differences can be found for structures with large unit cells, and approximately corresponds to the number of atoms in the system (Fig. 11). Considering all structures of the benchmark data set, both tools are similarly fast on average,

with LifeSoaks being 1.36 times faster than MAP_CHANNELS. For 1281 structures, MAP_CHANNELS is faster than LifeSoaks. However, the run-time difference exceeds 60 s in only 46 cases. The better average run time of LifeSoaks can be attributed to crystals with unit cells with larger numbers of atoms. When only considering systems with 50 000 and more atoms in the unit cell (12% of the data set), LifeSoaks is on average nearly 2.7 times faster than MAP_CHANNELS. This trend can be explained by the run time increasing roughly linearly [expected ] with the number of atoms with LifeSoaks, while it increases quadratically with MAP_CHANNELS. The latter has already been discussed in Section 1, outlining this major limitation of grid-based approaches.

Figure 11 Comparison of the run times of LifeSoaks and MAP_CHANNELS. (a) Comparison of the run times of both tools. (b) Dependence of the run time on the number of atoms in the unit cell for LifeSoaks. (c) Dependence of the run time on the number of atoms in the unit cell for MAP_CHANNELS.

Whereas the MAP_CHANNELS approach can be beneficial for coarse-grained analyses of solvent channels to obtain a crude estimate of the bottleneck radii, using LifeSoaks is preferable over using MAP_CHANNELS for the application scenario of predicting the bottleneck radii of crystals with small molecules to assess their applicability to soaking-based screening studies. Especially for particularly large unit cells, a concise assessment of the bottleneck radii is much faster than an analysis with MAP_CHANNELS. The availability of LifeSoaks on the ProteinsPlus web server enables on-the-fly soaking predictions without the hurdles of installation and parametrization. In addition, it reports the positions of bottlenecks as Cartesian coordinates, facilitating the detection of soaking obstacles for users. Finally, LifeSoaks offers automated analyses of binding-site accessibility, a feature that is not implemented in MAP_CHANNELS.