Engineering enzyme access tunnels

Enzymes are efficient and specific catalysts for many essential reactions in biotechnological and pharmaceutical industries. Many times, the natural enzymes do not display the catalytic efficiency, stability or specificity required for these industrial processes. The current enzyme engineering methods offer solutions to this problem, but they mainly target the buried active site where the chemical reaction takes place. Despite being many times ignored, the tunnels and channels connecting the environment with the active site are equally important for the catalytic properties of enzymes. Changes in the enzymatic tunnels and channels affect enzyme activity, specificity, promiscuity, enantioselectivity and stability. This review provides an overview of the emerging field of enzyme access tunnel engineering with case studies describing design of all the aforementioned properties. The software tools for the analysis of geometry and function of the enzymatic tunnels and channels and for the rational design of tunnel modifications will also be discussed. The combination of new software tools and enzyme engineering strategies will provide enzymes with access tunnels and channels specifically tailored for individual industrial processes.


Introduction
The textbook lock-and-key model for enzyme catalysis was first introduced by the Nobel laureate organic chemist Emil Fischer in 1894 (Fischer, 1894). The induced-fit and the selected-fit models, sometimes also referred as conformational selection, expanded the Fischer's rigid model to cover the flexibility of both the ligand as well as the enzyme (Chakraborty and Di Cera, 2017;Johnson, 2008;Koshland, 1995). The ligand binding induces conformational changes in the enzyme in the induced-fit model, whereas in selected-fit models the enzyme exists in equilibrium between multiple conformational states, some of which the ligand can bind to and stabilize the reactive conformation. Despite the old base theories, both induced and selected fit models describe well the catalysis of enzymes with the active sites localized on the protein surface ( Fig. 1). However, the majority (> 60%) of enzymes have their active sites located in deep internal cavities connected to the bulk solvent by access pathways (Marques et al., 2016;Pravda et al., 2014). These pathways are represented structurally by channels and tunnels. Channel or pore means a pathway leading throughout the protein structure, without any interruption by an internal cavity, with both sides open to the surrounding solvent. Tunnel means a pathway connecting a protein surface with an internal cavity or a pathway connecting more than one internal cavity. Tunnels enable the transport of substrates, cofactors, solvents, and products to and from the active site, and they play a crucial role in all parts of the enzyme's catalytic cycle. To include the impact of these access pathways to the biocatalysis, we have proposed a keyhole-lock-key model which better describes the catalysis of the enzymes with buried active sites and access tunnels (Marques et al., 2017a;Prokop et al., 2012).
The keyhole-lock-key model incorporates the passage of the ligands through the tunnels (keyholes) to the catalytic site of the enzyme and their exit from the site to the surrounding environment (Prokop et al., 2012;Wheeldon et al., 2016). Fig. 2 shows a simplified and generalized view of the enzymatic cycle. The catalytic cycle starts by (i) the passage of the ligand through the tunnel (Fig. 2B), (ii) followed by the reorganization of waters (Fig. 2C) and (iii) the binding to the catalytic residues (Fig. 2D). We note that the water dynamics and its impact on biocatalysis can vary from the one presented in the figure (Fink and Syrén, 2017;Hendil-Forssell et al., 2015;Kim et al., 2018aKim et al., , 2018b. The physical steps are followed by one or more chemical steps (iv) during which the product(s) are formed (Fig. 2E). The exit of products follows the opposite order from the binding, comprising of the physical steps of (v) the products unbinding from the catalytic residues with active site solvation ( catalysis alter the hydrogen bonding networks and/or the interaction of the protein with the surrounding solvent (Fields et al., 2015). The dynamics-based regulation of enzymatic catalysis tunes the biological function without affecting the structure of the ground state and modifications to different dynamic parts of the enzyme can uncover a spatial separation of the control of key enzymatic parameters (Saavedra et al., 2018). A combination of experimental and computational studies are often needed to study the dynamic properties of enzymes (Bottaro and Lindorff-Larsen, 2018).
The tunnels and loops as dynamic gating elements and the effect of water molecules on protein function are starting to attract the interest of enzyme engineers (Fink and Syrén, 2017;Hendil-Forssell et al., 2015;Jenkins et al., 2018;Kreß et al., 2018). The published experimental and computational studies clearly demonstrate that the engineering of the enzyme access tunnels has proven applicable for tens of enzymes from all the six enzyme classes. In the following sections, we will describe the success stories of enzyme access tunnel engineering both by humans and by nature, with special emphasis in the advances of the last 5 years, regarding engineering: (i) the catalytic activity, (ii) engineering the substrate and reaction specificity, (iii) engineering the enantioselectivity, (iv) engineering stability, and (v) the tools that can be used in the engineering of enzyme access tunnels. The discussed enzymes are shown in Table 1.

Engineering catalytic activity
Excluding the chemical reaction step (Fig. 2), all stages of the catalytic cycle of an enzyme with a buried active site utilize the access tunnels in either the transport of the solvent, substrates, cofactors or products. Many times, the catalytic rate is limited by either substrate binding or product unbinding (Cleland, 1975;Kreß et al., 2018;Narayanan et al., 2018). In these cases, the catalytic activity can be modified by adjusting the throughput of the ligands through the access tunnels. Besides the geometry of the tunnels, also the electrochemical and dynamic properties of the tunnels can be modified to ease the transport of ligands. Illustrative examples of how the tunnels can be modified to affect the enzymatic activity are shown in Fig. 4.
Blocking the transport tunnels by bulky residues is a common way to lower the activity of enzymes. Nature provides us with an example of human carbamoyl phosphate synthetase 1, which, unlike most other carbamoyl phosphate synthetases, cannot use glutamine as its source of ammonia (de Cima et al., 2015). This is partially due to the two ammonia transfer tunnels being blocked either by a Gln residue or by charged residues around the entrance. By enzyme engineering, Farrugia et al. modified tunnel lining residues in an accessory protein UreD to block the transport of nickel to a urease (Farrugia et al., 2015). The variants with blocked tunnel showed diminished urease activity, thus confirming the authors' hypothesis that the tunnel of UreD is the main route for the nickel transport.
Another interesting example of a blocked transport tunnel is shown with indoleamine 2,3-dioxygenase in which a Ser167His mutation coincides with a water tunnel leading to the active site (Lewis-Ballester et al., 2018). The tunnel is blocked by a hydrogen bond made by the introduced histidine residue to the opposing side of the tunnel. This leads to about 5000-fold reduction in the O 2 binding constant. However, the effect can be reversed with 3-indole ethanol binding, which results in a movement of a Phe side chain and a reorganization of the active site residues to the conformation represented by the open water tunnel.
The enzymatic activity can sometimes be affected by adjusting the tunnels of the interacting proteins and not the enzyme itself, such as in the case of lecithin-cholesterol acyltransferase (Gorshkova et al., 2018). The two Arg residues are positioned at the edge of a hydrophobic tunnel in the interacting dimeric protein ApoA-I. They help to position the fatty acids from the hydrolysis for the following esterification to cholesterol reaction. In their experiment, Gorshkova et al. mutated these Args of the ApoA-I to Ala and Gln, and both mutations lead to a significant reduction of the catalytic efficiency of the lecithin-cholesterol acyltransferase.
Controversially, it is also possible to increase the catalytic activity by introducing bulky mutations into the tunnels. Pavlova et al. introduced bulky mutations into access tunnels of haloalkane dehalogenase DhaA by directed evolution and ended up with a 32-times more active mutant towards the anthropogenic substrate 1,2,3-trichloropropane (Pavlova et al., 2009). They rationalized that the increased activity is due to the occluded active site promoting the active conformation of the substrate because fewer water molecules are diffusing to the active site and the substrate has a more limited conformational space in the active site. Later, molecular dynamics simulations showed that high catalytic efficiency is due to the occluded site together with a decrease in the energy barrier during the reaction step (Marques et al., 2017b).
Logically, the enzymatic activity can be also affected by replacing bulky residues with smaller ones, which might result in better transport of substrates, products and modified solvation of the active site. By mutating bulky residues Leu, Met and Phe in the access tunnels to smaller Ala Kong et al. managed to increase the enzymatic activity of an epoxide hydrolase up to 42 times (Kong et al., 2014). The rate of the product release was increased significantly and also the substrate specificity was shifted towards the bulky epoxides.
Brezovsky et al. combined the effects of both approaches, firstly, closing the main access tunnel, and secondly, opening an auxiliary access tunnel by modifying a haloalkane dehalogenase LinB . This mutant proved to be the most catalytically efficient haloalkane dehalogenase isolated or constructed in the laboratory. The mutant with extensively re-engineered access pathways showed a broad substrate specificity. The resulting high enzymatic activity is at least partially caused by the increased ligand transport due to the de novo tunnel, and it is further facilitated by the changed dynamics of the main tunnel .  3. Examples of access tunnels detectable in the crystal structures of enzymes cytochrome P450 and carbamoyl phosphate synthetase. A) Multiple access tunnels (various colours) leading to the active site of the human cytochrome P450 CYP3A4 (PDB 5VCC). B) A long tunnel (yellow) connecting two active sites (magenta) in human carbamoyl phosphate synthetase (PDB 1A9X). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1
An overview of studies where the tunnels of enzymes have been modified using the protein engineering methods. These articles describe the same enzyme but to preserve the nomenclature presented in the references they are shown as the separate studies.
P. Kokkonen, et al. Biotechnology Advances xxx (xxxx) xxx-xxx The dynamic nature of tunnels is impacting their function. This is clearly seen in a glycoside hydrolase, chitinase B, which forms a tunnellike catalytic cleft in which the chitin polymer chain can move upon substrate binding (Hamre et al., 2017). This enables the formation of new reactive complexes and a processive reaction chain. Four mutations of the tunnel-lining residues increased the conformational freedom of the tunnel walls, implying that these residues likely play the key roles in the structural rigidity and the shaping of the tunnel before the substrate binding. Hamre et al. displayed that the tunnel properties are crucial for both the substrate binding and the processivity of the enzyme (Hamre et al., 2017).
The dynamic tunnels can also stabilize the active conformation, as is the case of a histone deacetylase 3. The binding of a cofactor to this enzyme restricts its conformational space to a rigid conformation which primes the enzyme for its substrate recognition (Arrar et al., 2013b). A histone deacetylase 3 variant lost all of its deacetylase activity towards the natural substrate due to a single Arg to Pro mutation in its access tunnel, even though the variant adopted a conformation similar to the rigid ready-to-bind substrate conformation in the wild-type (Watson et al., 2012). Though the variant stabilizes the rigid and activated conformation, the authors show through molecular dynamics simulations that the active conformation is not as stable as in the wild-type because the tunnel is fluctuating (Arrar et al., 2013a;Watson et al., 2012). In another study, Medlock et al. constructed an inactive variant of a ferrochelatase by modifying a residue from a Phe to Arg in a junction of two access tunnels and the active site (Medlock et al., 2012). This disrupted the hydrogen bond network and the reoriented the structurally important π helix, thus interfering with the product release. Fig. 4. Schematic representation of the modifications of access tunnels for engineering enzymatic activity. (A) Replacing small residues with bulky residues to block a tunnel and hinder ligands or solvent transport. (B) Replacing bulky residues with smaller residues to make the tunnel larger or (C) to create a new tunnel, increasing the ligands or solvent transport. (D) Affecting the dynamics of the tunnel by modifications of dynamic elements. Enlarging the tunnel generally increases enzyme promiscuity, letting also larger ligands pass. (C) The shape of the tunnel can be changed to let specific ligands pass. P. Kokkonen, et al. Biotechnology Advances xxx (xxxx) xxx-xxx

Engineering substrate specificity and promiscuity
Enzymes have evolved to catalyze one or more specific reactions, which is one of their benefits over the far less specific chemical catalysts. The buried active sites are only accessible by substrates which can pass through the enzyme access tunnels. Thus, these tunnels offer an excellent engineering target for differentiation of substrates (Fig. 5). In addition to controlling a substrate access, the tunnels also contribute to the reaction specificity in case of multifunctional enzymes.

Substrate specificity
One of the widely applicable methods to change the substrate specificity is to shorten or lengthen the tunnel that binds the carbon tail of fatty acids or similar long-chained compounds in enzymes that use these molecules as substrates. For example, a fatty acid dehydratase FabA has a tunnel which binds the acyl chain of the substrate, and the length of this tunnel determines the substrate specificity of this enzyme (Finzel et al., 2015). The aldehyde-deformylating oxygenase, applicable in the biofuel production displays low specificity among substrates with various carbon chain length (Bao et al., 2016). Bao et al. introduced bulky modifications narrowing the access tunnel and this resulted in a preference for the shorter chain aldehydes, making the enzyme more lucrative target of the biofuel industry. Similarly, by modifying a long tunnel with a bulky residue, an alkane hydroxylase could be made to prefer substrates with a shorter chain (van Beilen et al., 2005). Mutations blocking a part of the tunnel seems to be a general strategy to allow binding of the ligands of particular lengths.
Homology modelling was employed to design Asp131Ser mutation in an access tunnel of a cytochrome P450 CYP153A35 (Jung et al., 2018). This resulted in a 17-fold improvement in the total turnover number and catalytic efficiency towards palmitic acid hydroxylation in vitro, but only 50% improvement of palmitic acid yield in the whole cells. With computational docking experiments, the authors demonstrate that Ser residue forms hydrogen bonds most effectively with palmitic acid-sized compounds. The resulting effect of increased in activity is synergistic with mutations in the active site.
The blockage of an access tunnel can lead to a viral relapse in patients suffering from a hepatitis C infection with a mutation introduced to RNA-dependent RNA polymerase NS5B (Uchida et al., 2018). Here, the mutation Ala218Ser causes the tunnel to hinder the binding of the inhibitor sofosbuvir, which prevents its conversion into the active metabolite, thus leading to inefficient treatment. The mutation does not significantly affect the natural function of the enzyme, as the passage of the nucleotide triphosphate is not affected as much as the inhibitor.
A membrane-bound fatty aldehyde dehydrogenase employs an intricate mechanism to control the entrance of the substrates to the access tunnel (Keller et al., 2014). It has a helix which can move to block the access tunnel and stabilize the ligand in the reactive conformation. A mutation in this helix region causes Sjögren-Larsson syndrome, an autosomal neurocutaneous disease. Interestingly, the removal of this gating helix did not compromise the overall catalytic activity but shifted substrate specificity towards compounds with shorter chains.
Substrate binding to α-ketoglutarate dependent oxygenase factor inhibiting HIF limits the solvent access to the active site (Chaplin et al., 2018). By expanding this tunnel with Ala replacements, the substrate inhibition by α-ketoglutarate increased. Since the mutations did not affect the binding of molecular oxygen to the active site as the authors expected, they rationalize that the narrowing of the tunnel caused by the substrate binding prevents the larger molecules from reaching the cofactor once the first substrate has bound. Subramanian et al. conducted an extensive screening (> 50 variants) of pig and human D-amino acid oxidase water access tunnels (Subramanian et al., 2018). By facilitating the solvent and substrate access they affect both the activity as well as the specificity of engineered enzymes. By identifying the tunnels and analyzing the water passage through tunnels from molecular dynamics simulations, they could pinpoint the most significant changes to residues Tyr55, Tyr224 and Tyr314. These residues are located on different loops around the entrance of the access tunnels. Especially the mutations of Tyr55 contributes to substrate specificity and the mutation Tyr55Ala increased the activity towards bulky substrates. This residue is important for anchoring the loop that covers the tunnel to active position and separating the exits of the two analyzed tunnels.

Reaction promiscuity
Many enzymes catalyze multiple chemical reactions and their rates are affected by the enzyme tunnels controlling the access of substrates, solvent or cofactors to the active site. By modifying a residue in the active site and three residues at a water tunnel, David et al. made an agarase from Zobellia galactanivorans lose most of its hydrolase activity while retaining and enhancing its trans-glycosylase activity . This increased the water flux into the active site, which counterintuitively reduced the hydrolase activity. The authors hypothesize that the hydrolytic water molecule needs to form a favorable network of hydrogen bonding interactions to anchor it in a favorable geometry for the hydrolysis, and this is hindered by the increased water dynamics in the active site.
Another example of a promiscuous enzyme is hydrolase Mhg which displays both γ-lactamase and perhydrolase activities despite having its identical catalytic triad and active site geometry to that of an esterase (Yan et al., 2016). Strikingly, by mutating Leu in the substrate binding tunnel into smaller residues the researchers removed the steric blockage of the tunnel that prevented the entrance of ester compounds and thus introduced the esterase activity. The designed esterase variants catalyzed in addition to classical substrates a wide range of diverse chiral substrates, even though there were differences in their reaction preference. Despite starting from a promiscuous enzyme, the variants from the saturation mutagenesis retained only a single kind of activity: (i) the Leu233Met variant had perhydrolase activity, (ii) the Leu233Thr variant had γ-lactamase activity and (iii) the Leu233Gly variant had esterase activity. Clearly, this site provides a functional hot-spot for the enzymatic evolution of different activities, and there is a trade-off between the different catalytic activities as not a single variant retained all the three activities. Syrén et al. used an ingenious approach to engineer a more catalytically effective triterpene cyclase from Alicyclobacillus acidocaldarius (Syrén et al., 2014). They identified three water channels leading to the active site, and knowing that the substrate binding is hindered by a high entropic cost of the pre-folding of the substrate for cyclization, they wanted to increase entropy by releasing bound water molecules from these water channels. By blocking the water channels with Trp point mutations (originally Ser or Phe) they could get up to 100-fold higher apparent k cat /K M values for a pentacycle formation. They also observed close to perfect enthalpy-entropy compensation, indicating that different extents of cyclization could operate by the same reaction mechanism. In the continuation study, the molecular dynamics simulations of the tunnel blocked mutants unexpectedly revealed increased water flow around the catalytic Asp residue (Gustafsson et al., 2017). The authors hypothesize that the increased water flow switches the reaction mechanism towards a Grotthuss-like proton transfer, which is characteristics by an increased protonation distance and negative activation entropy.

Engineering enantioselectivity
The production of pure enantiomers or the utilization of single enantiomers as substrates is a highly valuable property of enzymes. Lipases are a family of enzymes which contain a tunnel region that binds the acyl end of their substrates (Kingsley and Lill, 2015). The length (from 8 to 22 Å) and geometry of this tunnel vary a lot among different lipases (Pleiss et al., 1998). By engineering the tunnel-lining residues one can affect the specificity for the length of the acyl chain (Bassegoda et al., 2012;Juhl et al., 2010;Schmitt et al., 2002), as discussed in the substrate specificity section. Moreover, the geometric properties of the tunnel can be modified to differentiate saturated and trans fatty acids from cis fatty acids because the former can adopt linear geometries whereas the latter cannot (Fig. 5C) (Brundiek et al., 2012). Selectivity towards saturated and trans fatty acids was accomplished by the introduction of a narrow region with bulky residue replacements. Conversely, cis-selective lipases have tight curvatures in their access tunnels to prevent the entrance of the long and straight trans fatty acids.
Li et al. combined the strategy of tunnel and active site modifications to develop monoamine-oxidase MAO-N variants with different enantioselectivities . The variants that combined both tunnel-lining and active site mutations proved to be more enantioselective and active towards larger and more hydrophobic substrates. The directed evolution method decreased the polarity of the access tunnel, making the entry and egress of hydrophilic substrates easier. The mutations targeting active site and the tunnel-lining residues were needed since the variant carrying a single tunnel mutation was found inactive.
The oxygenation position and stereospecificity of lipoxygenases can be altered by an Ala to a Gly mutation of a residue lining the active site and the oxygen access channel (Coffa and Brash, 2004;Kalms et al., 2017;Suardíaz et al., 2013). This is seen in the specificities of the lipoxygenases from different organisms: the S-lipoxygenases have a conserved Ala residue in this position, whereas R-lipoxygenases have the smaller Gly. The molecular mechanisms of this change in hydroperoxidation reaction of arachidonic acid were thoroughly explored by Saura et al. (Saura et al., 2017). In the wild-type lipoxygenase of that study, the reaction occurs on the C 8 of arachidonic acid with an R stereochemistry. The Ala mutant blocks the oxygen access tunnel in the region close to C 8 and redirects the oxygen towards C 12 , where the oxygenation happens with an S stereochemistry.
Liskova et al. demonstrated in a study of two haloalkane dehalogenases that the similar enantioselectivity of the enzymes was caused by different factors (Liskova et al., 2017). Whereas the naturally highly enantioselective DbjA has an open, solvent accessible active site, the engineered enzyme DhaA31 has a narrow tunnel connecting the active site with the solvent. DhaA31 was engineered from low enantioselectivity DhaA with the intention to increase enantioselectivity. In the end, the DhaA31 and DbjA display similar enantioselectivity when converting 2-bromopentane. The enantioselectivity in the DhaA31 variant is caused by two out of five tunnel and active site substitutions, which caused the preferred substrate to bind in a favored configuration. On DbjA, the enantioselectivity is explained by preferred binding mode due to water-promoted interactions.

Engineering protein stability
Thermodynamic and kinetic stabilities are crucial properties of enzymes to be optimized for industrial applications (Iyer and Ananthanarayan, 2008). According to the widely accepted Lumry-Eyring model, the inactivation of enzymes consists of at least two steps: (i) the reversible unfolding of the native state enzyme and (ii) the kinetically irreversible steps leading to changes in the crucial stabilizing interactions and/or to enzyme aggregation (Andrews and Roberts, 2007;Li and Roberts, 2009). Many factors, such as temperature, pH, proteases, solvents, other chemicals and ionic strength affect the stability of enzymes. Access tunnels are formed in such a way that ligands and solvents can pass to the catalytic site, they create voids and decrease a number of intra-molecular interactions within the protein structure, which generally decreases the conformational stability. Thus, reducing the size and/or the number of the access tunnels generally increases the stability of enzymes (Burley and Petsko, 1985;Georis et al., 2000;Meurisse et al., 2003;Serrano et al., 1991). However, this strategy to increase the protein stability is not so commonly used. The lack of studies might be because this engineering strategy is naturally prone to stability-activity trade-off and a careful optimization of mutations inside the tunnels is needed to enhance stability without scarifying catalytic activity.
By removing 4 or 10 residues of the N-terminal loop structure, which gates the access of substrates to the active site, Singh et al. lowered the thermal stability of a hyperthermostable carbohydrate esterase by 30 to 40°C of the optimal reaction temperature (Singh et al., 2017). The modified N-terminal loop is assumed to control the specificity of the substrate entrance to a 10-Å long access tunnel (Vincent et al., 2003). However, only a modest preference for acetylated xylose over acetylated glucose was observed. The authors concluded that the N-terminal loop stabilizes a rigid near-by loop that harbors one of the catalytic residues, thus it was selected by evolution to increase the stability, in addition to its modest role in substrate selectivity. Koudelakova et al. mutated four residues in the access tunnels of a haloalkane dehalogenase which resulted in a 19°C increase in the melting temperature (Koudelakova et al., 2013). The change also extended the half-life of the enzyme in 40% dimethyl sulfoxide solution from minutes to weeks. With protein crystallography and molecular dynamics simulations, the authors concluded that the increases in protein stability are caused by more tight packing of one of the protein domains, i.e. more constricted access tunnels. The four mutations, however, caused a decrease in enzymatic activity. In a follow-up study, the researchers improved the activity by mutating forth two of the four previously modified tunnel-lining residues (Liskova et al., 2015). Finally, only one mutation in the access tunnel was shown to retain almost all of the thermostability (4°C difference) of the 4-point mutant while simultaneously showing higher enzymatic activity with the substrate specificity of the wild type.
The stability of lipases is limited in many industrial processes by the denaturation of polar alcohols present in the reaction mixture. To increase the stability in methanol of a lipase from Geobacillus stearothermophilus, Gihaz et al. introduced new aromatic interactions into the long tunnels leading to the enzyme active site (Gihaz et al., 2018). They were careful to maintain native hydrogen bonding networks and changes to the catalytic and metal binding residues. The three stabilizing mutations found in the initial screening were exclusively in the longest access tunnel leading from the solvent to the active site. The continuation of the study with doublepoint combinations of the original mutations revealed that the relationship between stability in methanol, melting temperature (T m ) and activity is rather complex. The most increased stability, T m value and activity were displayed by the double-point mutant Ala187Phe + Leu360Phe, whereas the second double-point mutant Ala187Phe + Leu184Phe displayed increased stability and T m , but lower activity. The third double-point Leu184Phe + Ala187Phe displayed decreased methanol stability, lower T m and lowest catalytic activity, probably due to significant aromatic cluster rearrangement, which disturbed the overall stability of the protein.

Software tools for tunnel analysis and engineering
As the software tools for the determination of enzymatic tunnels have been reviewed by Brezovsky et al. in 2013 andLill in 2015, and the tools for protein engineering have been reviewed by Ebert and Pelletier in 2017, only a brief overview of the current state of the software tools is given here Ebert and Pelletier, 2017;Kingsley and Lill, 2015). The access tunnels and channels can be identified with geometric methods, such as CAVER, MolAxis, MOLE, CHEXVIS, 3 V and CCCPP (Table 2) (Abrams and Bussi, 2013;Chovancova et al., 2012;Masood et al., 2015;Voss and Gerstein, 2010;Yaffe et al., 2008). The downside of most of these methods is that they are designed to find tunnels in individual static structures, and generally the dynamic nature of tunnels is forgotten (Chovancova et al., 2012;Kingsley and Lill, 2014a).
The dynamic nature of protein tunnels can be studied with NMR experiments or molecular dynamics simulations (Boehr, 2018;Campbell et al., 2018;Kleckner and Foster, 2011;Zhao et al., 2018). Both of these methods have their limits: it is challenging to study large (> 50 kDa) proteins with NMR techniques, and the molecular dynamics simulations typically only cover timescales up to milliseconds whereas the dynamic movements of proteins and their tunnels might occur in timescales of seconds or longer (Boehr, 2018). To combat the problem of low timescales sampled in molecular dynamics studies, researchers conduct simulations in elevated temperatures or use enhanced sampling techniques such as accelerated molecular dynamics or metadynamics (Abrams and Bussi, 2013;Salmaso and Moro, 2018). The resulting huge ensembles of protein and tunnel conformations from NMR studies or simulations require a form of clustering to separate the tunnels from each other, such as the one used in CAVER 3.0 (Chovancova et al., 2012). The conformational ensembles of protein tunnels obtained from NMR studies or molecular dynamics simulations offer a possibility to find occasionally closed tunnels and to discover various tunnel geometries ( Fig. 6) (Chovancova et al., 2012;Jurcik et al., 2018). Caver Analyst 2.0 can be used to study the opening of an access tunnel of a dehalogenase during a molecular dynamics simulation as demonstrated in Fig. 6B and C.
The software tools described above analyze primarily the geometrical properties of tunnels with possibly mapped physico-chemical properties, while they do not consider tunnels' functionality. The functionality of tunnels can be studied with random expulsion molecular dynamics simulations or IterTunel which uses steered molecular dynamics together with a geometric tunnel calculation algorithm (Lüdemann et al., 2000;Kingsley and Lill, 2014b). These molecular dynamics methods are computationally demanding which is why there are a few software tools to study ligand binding that use less demanding molecular docking or robotic algorithms. Moma-LigPath, SLITHER, ART-RRT and CaverDock are the software tools based on molecular docking and/or robotics suitable for analysis of the passage of ligands in and out of buried active sites (Table 3). One of the most recent tools, CaverDock implements a novel algorithm which is based on molecular docking and can produce contiguous ligand trajectory and estimation of the binding energy along the pathway (Fig. 7) (Filipovič et al., 2019). CaverDock, and the other software tools in Table 3, can be used in comparing the tunnel function between different variants of an enzyme, as well as comparing the throughput of ligands between different tunnels in the same enzyme. All of these tools work only on static protein structures.
There is also a database which contains information about enzymatic tunnels: tunnels calculated for X-ray crystal structures can be found in ChannelsDB database, which provides information about the position, shape, and properties of tunnels and channels of the proteins covered by the Protein Data Bank (Rose et al., 2013;Pravda et al., 2018a). The database offers 3D visualization of the enzyme of interest and its tunnels, as well as annotations for individual residues retrieved from the UniProt database ("UniProt,", 2017).
Once the tunnel-lining "hot spots" have been identified by geometrical approaches, by the analysis of ligand passage or taken from the databases, the engineering of tunnels continues with the prediction of the effects of the mutations to the enzyme structure. Recently, a detailed computational protocol employing CAVER was published to simplify this process; the analysis of tunnel lining residues and the selection of "hot spot" residues for tunnel engineering Chovancova et al., 2012). There are many web servers and standalone software tools for predicting the effect of point mutations on protein function and stability (Table 4). HotSpotWizard 3.0 helps the users to identify the functional and stability "hot spots" and to design libraries or mutations for experimental characterization (Sumbalova et al., 2018). It is also part of the newly launched collection of methods for in silico protein engineering available via NewProt portal (Schwarte et al., 2017). The portal offers access to many enzyme engineering tools for the introduction of single-point and multiple-point mutations. The portal uses HOPE software to indicate mutations with possible undesirable side effects, such as loss of activity or stability (Venselaar et al., 2010). Another software, the computer-aided directed evolution of enzymes (CADEE), which prepares and analyzes in silico semi-automated directed evolution of enzymes. CADEE 0.9.1 reduces the screening efforts and tries to maximize the likelihood of identification of successful mutants (Amrein et al., 2017). Selected mutations are analyzed by an empirical valence-bond approach which combines the computational efficiency of the force-field descriptors with a quantummechanical description of chemical reactivity (Warshel and Weiss, 1981). The software runs the computational mutagenesis of the userdefined residues using empirical valence-bond simulations, thus reducing the experimental work because only the most promising mutants need to be experimentally characterized. Other web services provide help when designing stable mutants with modifications in the access tunnels or channels. FireProt web server offers users to design thermostable multiple-point enzyme variants with a combined strategy of evolutionary and energetic methods Musil et al., 2017). Protein Repair One Stop Shop (PROSS) web server offers an algorithm to design stable proteins that can be expressed in E. coli (Goldenzweig et al., 2016).

Conclusions and perspectives
The case studies described in this review demonstrate that the modifications to the enzyme access tunnels and channels provide a feasible expansion to the active site engineering strategy. By modifying residues in the access tunnels and the active site, it is possible to affect practically all the stages of the catalytic cycle of the enzyme. The optimization of the access tunnels is an emerging protein engineering strategy with immense possibilities, which are seen in the steep rise in the number of studies regarding the area in the last couple of years. The field is still developing as further insights into the physico-chemical properties and dynamics of the access tunnels are being revealed. These insights provide the protein engineers with new sites for the modifications, leading to changes in important biocatalysts' properties: (i) catalytic activity, (ii) substrate specificity, (iii) enantioselectivity and (iv) protein stability.
Various computational methods can be used for the detection and analysis of enzyme tunnels and for identification of tunnel-lining, bottleneck or gating residues. Predictions of specific mutations, which will lead to desired changes in catalytic activities, substrate ranges or stabilities are still challenging. Computational design decreases the amount of experimental work in protein engineering which should ultimately lead to the rational design of the enzyme access tunnels suitable for specific ligands. Development of the software tools for identification of tunnels in proteins and other biomolecules has received momentum during the last decade during which a large number of tools have been introduced to the community. Majority of the tools identify tunnels in the static structures. Caver 3.0 and Caver Analyst 2.0 enable the analysis of dynamic structures and can handle even very large ensembles of structures from molecular dynamics simulations. A computational analysis of these ensembles is challenging but essential since Fig. 6. Analysis of protein tunnels from molecular dynamics simulations using Caver Analyst 2.0 (Jurcik et al., 2018). A) The two tunnels identified in a single snapshot in a haloalkane dehalogenase during a 50 ns (500 snapshots) molecular dynamics trajectory. B) The cross-section of the bottleneck in the blue tunnel during the simulation and the surrounding residues (coloring: dark -the start of a simulation, orange -the end of a simulation). Each line represents the cross-section of a single snapshot. C) Tunnel profile of the blue tunnel during the simulation (coloring: dark -the start of a simulation, orange -the end of a simulation). Each line represents the tunnel profile in a single snapshot. The vertical line indicates the cross-section displayed in B. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 3
Software tools for studying the transport of ligands in the access tunnels. Another advancement in the field is the computational analysis of ligand passage through protein tunnels. This is an essential complement of geometrical analysis, which precisely identifies a location of access/ egress pathways within protein structures, but do not provide information whether these pathways can be used by cognate ligands. Continuously growing computer power and enhanced sampling methods of molecular dynamics allow an analysis of ligand binding and unbinding at the atomistic details. New tools like Moma-LigPath or CaverDock provide information on ligand binding trajectories. These applications are yet to be developed by the community of structural biologists, enzymologists and medicinal chemists. Further improvement of existing algorithms to include for example protein dynamics and achieve higher preciseness in estimating energy barriers without comprising the speed of calculation. Development of these tools is essential for achieving the ultimate goal of designing protein tunnels optimized for specific cognate ligands. Fig. 7. Analysis of ligand transport through a tunnel using CaverDock. A) Dihydroorotic acid passing through the access tunnel of dihydroorotate dehydrogenase (PDB ID 1D3G; Liu et al., 2000). B) The energetic profile of the ligand binding and the tunnel profile. The colored dots along the energetic profile indicate the corresponding snapshots displayed in A. The ligand transport analysis displays a clear energy barrier in the access pathway around the red ligand. The tunnel geometry displays similar bottlenecks also closer to the surface, even though they do not display such large energetic barriers. This can be explained by the oval shape of the tunnel at those points: the radius is determined by the largest ball that can fit through, whereas a cyclic ligand, such as dihydroorotic acid, can adjust its orientation to pass. CaverDock requires tunnels calculated by Caver and approximates the description of the chosen tunnel with a set of discs. CaverDock then docks a ligand along the tunnel, constraining one atom to a disc at a time. The binding energies, obtained using the AutoDock Vina algorithm, are then plotted for every step along the tunnel profile (Filipovič et al., 2019). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 4
Software tools to assist in engineering protein access tunnels. Design of stable and soluble proteins (Goldenzweig et al., 2016) P. Kokkonen, et al. Biotechnology Advances xxx (xxxx) xxx-xxx