{'fraction_non_alphanumeric': 0.034146652226540106, 'fraction_numerical': 0.013505765432885265, 'mean_word_length': 4.642703091301222, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A new method for finding the minimum free energy pathway of ions and small molecule transportation through protein based on 3D-RISM theory and the string method

A new method for finding the minimum free energy pathway (MFEP) of ions and small molecule transportation through a protein based on the three-dimensional reference interaction site model (3D-RISM) theory combined with the string method has been proposed. The 3D-RISM theory produces the distribution function, or the potential of mean force (PMF), for transporting substances around and inside proteins. By applying the string method to the PMF surface, one can readily determine the MFEP of the transporting substances inside a protein. The method has been applied to consider the Na+ conduction pathway of channelrhodopsin as an example.


Introduction
The transportation of ions and molecules through proteins is one of the most fundamental life phenomena to maintain living systems [1,2]. An ion channel adjusts the ion concentration according to the difference of the chemical potential of both intra-and extracellular sides of the cell membrane. It is also related to the signal transduction processes and the conduction of electrical signals. While channels allow the transporting of ions or small molecules across the membrane only passively, a pump protein transports the substances actively across the membrane against chemical potential gradients. Because of active transport, pump proteins create membrane potentials or proton gradients between the intra-and extracellular sides of membranes. [3] The characteristic of ion (or small molecule) transportation is expressed by the free energy profile or potential of mean force (PMF) through the protein. For example, the barrier on the PMF is related to the ion conduction flux. Therefore, the evaluation of the PMF is essential to understanding the mechanism of the transportation. However, computation of the PMF is a challenging issue for computational science because of the difficulties with evaluation of the PMF acting on the ions.
The PMF is defined as the potential that gives the average force taken for configurations of protein, water, and ions. Because ion transportation is highly coupled with the structural fluctuation of the protein and membrane, and the distribution of solvent species, the configurational integral for the PMF should include a sufficient number of configurations of both the protein structures and solvent distributions. Therefore, efficient sampling methods of the protein structures and the distributions of solvent, including water and ions, for the configuration integral are highly desirable. Many theoretical and computational approaches have been proposed to evaluate the PMF based on molecular dynamics (MD) simulations [4][5][6]. A nonphysical sampling such as umbrella sampling, metadynamics, and adaptive biasing force method is often employed to accelerate the MD sampling through the reaction coordinate of ion transportation [7][8][9]. However, these methods become inefficient when the number of reaction coordinates increases. Therefore, a method to determine the appropriate reaction coordinate is also important.
The difficulties in the sampling of the solvent distribution can be overcome by the threedimensional reference interaction site model (3D-RISM) theory. The spatial distribution function (SDF) of the solvent by 3D-RISM is obtained through a complete ensemble average over the entire configuration space of solvent molecules, including ions in the thermodynamic limits [10][11][12][13][14][15]. Therefore, by applying 3D-RISM, one can determine the PMF surface of ions for the entire space prior to the reaction coordinate search for a given protein structure. To find the reaction coordinates on a given potential surface, the string method is effective [16][17][18] and has been applied successfully to various problems.
In this paper, we propose a new method for finding the minimum free energy pathway (MFEP) of ions (or small molecules) on the PMF surface obtained by 3D-RISM theory combined with the string method. This method consists of three steps. First, the SDF of ions is evaluated by 3D-RISM theory for given protein structures. The protein structures are assumed to be taken from the MD trajectory. (However, we employ the protein structure taken from the protein data bank (PDB) for simplicity in the present paper.) The PMF is readily obtained from the SDF. Second, the stationary position of the ions inside a channel is determined based on the SDF, where the SDF shows local maximum values. Third, the string method is used for the search of the MFEP that connects two stationary points on the PMF surface determined in the previous step [19][20][21].
A similar method based on the 3D-RISM theory and the elastic band method was proposed by Tanwar, Sindhikara, and co-workers, which is referred to as RismPath [22]. Because the method proposed in this paper is an extension of RismPath, we refer to the new method as RismPath/SM.
The advantages of the string method over the elastic band method are higher stability and accuracy, which have been well discussed by E et al. [19]. The string method is a powerful tool to find an MFEP on a given surface that has been applied to various problems of finding the MFEP for a rare event involving a protein structure change [19,[23][24][25]. In contrast to those approaches, in the present method, the string method is used only for optimum route search that connects two points on the PMF surface with a given protein structure. In this study, the Na + conduction pathway of channelrhodopsin (ChR) is examined as an example to demonstrate the RismPath/SM method. In addition, to consider the effects of the protonation state of the amino acid residue on the pathway, both the protonated and unprotonated states of E136 and E140 have been examined. The MFEPs and the corresponding free energy profiles will be discussed.

Method
The computational procedure to determine the MFEP of the transporting substances through the protein is schematically depicted in Scheme 1. Although the present method can be applied to the MFEP search, not only for ions but also for small molecules, we focus on the ion case in this paper for simplicity.
First, to obtain the distribution function of the ions the 3D-RISM calculations should be conducted. Here, the target protein is regarded as a "solute" while the ions and water are treated as "solvents." The distribution function of solvent atom g, � ( ) is obtained by iteratively solving the 3D-RISM equation [10][11][12] and the Kovalenko-Hirata closure equation [13,26] where * denotes the convolution integral, ℎ � ( ) = � ( ) − 1 and � ( ) are total and direct correlation functions, and �� ( ) is a bulk solvent susceptibility, which is obtained by solving the RISM for the bulk solvent system prior to the 3D-RISM calculation. = 1/ B is the thermodynamic beta. The solute-solvent interaction potential � ( ) is given by where �� is the distance between solute atom a and solvent atom g. �� , �� , and � are the Lennard-Jones parameters and partial charge, respectively, with the usual meaning. Because the distribution functions obtained by 3D-RISM theory are obtained through a complete ensemble average over the entire configuration space of solvent molecules unlike MD, the distributions of ions inside a protein cavity or channel can be reproduced automatically. The PMF of solvent atom g is readily obtained by It is noted that the PMF acting on the transporting species should include not only the solvent distribution also the distribution of the force from the fluctuating protein. However, in the present method, since we assumed the fixed protein structure, the PMF can be expressed as Eq. (4).
Because the distribution function, � ( ), represents the probability of finding the solvent atom g at position , the probable position to accommodate the ions can be specified based on the distribution functions. The Placevent algorithm is a simple algorithm explicitly to place solvent molecules using the 3D-RISM distribution function [27]. To place the solvent molecules, the population function is defined as where � is a number density of solvent atom g in bulk solvent. The superscript of the population function denotes the number of iterations. The first explicit solvent atom is placed at the location with the highest population. After placing one explicit solvent atom, a new population function is obtained to maintain the total population of solvent atoms as where (�) denotes the area inside the sphere of radius � centered at � where the i-th explicit solvent atom is placed. � is determined to satisfy the condition: once the explicit atom has been placed and the corresponding distribution has been reduced. The next explicit atom will be placed based on the new population function. These processes are repeated until a satisfactory number of explicit solvent atoms are placed.
The position of the explicit atom placed above may correspond to a stable or metastable position of the transporting substrate atoms. In the present study, the MFEP connecting two stationary points of ions is determined by the string method [19,23,24]. This is a distinguishing point of the new method proposed here from the RismPath method; therefore, we call the new method RismPath/String Method (RismPath/SM) hereafter. In the present study, we employed the zero-temperature string method for simplicity [23,24]. It is noted that the finite temperature string method is also available, and it may be more effective in the case of a complicated PMF surface [20,21]. In the string method, a string connecting initial and end points of the path is discretized using N discretization points, which are referred to as images.

Computational details
To demonstrate the RismPath/SM method, we consider the Na + conduction pathway of ChR as an example. The protein structure is taken from the PDB entry code 3UG9, which is a crystal structure of the closed state of ChR (see Fig. 1A) [28]. Although the structure is a closed state, the pathway connecting the retinal-binding pocket to the extracellular region has been suggested, which is colored blue in Fig. 1B. In the present study, we used a fixed protein structure for simplicity. To consider the effects of the protonation state of the amino acid residue in the middle of the pathway, both the protonated and unprotonated states of E136 and E140 have been examined (see Table 1). We employed the CHARMM36 force field parameters for a standard amino acid [29] and the parameters for the Schiff base (SB) are taken from ref. [30]. The 3D-RISM calculations have been conducted using our own 3D-RISM program package [31]. NaCl solution (1 M) in ambient condition is assumed with the TIP3P water model at 298.15 K. The number of grid points in the 3D-RISM calculations was 256 3 with a spacing of 0.5 Å. For the string method, we used a modified version of the string method program originally coded by Ren [19,23]. Twenty-four images were prepared for each string, which corresponds to placing one image approximately every 1 Å. We employed Hermite interpolation to the PMF surface by 3D-RISM to obtain the PMF value at arbitrary coordinates for the string method.

Results and discussion
First, the 3D-RISM calculations for evaluating the distribution function of Na + , Cl -, and water have been conducted. The distribution functions of Na + and the oxygen of water around the ChR model (a) are depicted in Fig 2A and 2B, respectively. The isosurfaces of Na + ( ) = 1.2 and O ( ) = 1.2 are colored green and red, respectively. The noteworthy distribution of Na + can be observed around E136 and E140, whereas the distribution of oxygen of water spreads throughout the channel. A small Na + distribution can also be found near SB.
Based on the distribution functions, the explicit water and Na + are placed by the Placevent algorithm. In Fig 2C, some of the explicit atoms are depicted, which will be employed as the start, intermediate, and end points of the string method. The Na + atom labeled "S" is coordinated to the side chain of D292, and that labeled "I" is coordinated to the side chain of E136. "E A " and "E B " are located near the exit on the extracellular side. The oxygen atoms of water labeled "I A " and "I B " are located in I and E A , I and E B , respectively. respectively. (C) The explicit Na + , green, and water oxygen, red, placed using Placevent.
Here, we assumed two routes of the Na + conduction pathway, namely, S-E A and S-E B , which are referred to as routes A and B, respectively. To prepare the initial guess of the images on the string connecting S and E A or E B , the atoms "I," "I A ," and "I B " are used as the intermediate sites.  shown. The small spheres colored green and red denote the images on the MFEPs of routes A and B, respectively. Each image is placed at 1.09 Å and 1.11 Å for routes A and B, respectively.
The PMF profiles along the MFEPs are plotted against the images in Fig. 4, where the image 0 corresponds to the atom "S" while the image 25 is "E A " (or "E B "). At a glance, a peak can be found at images 2 to 3. This is because the positively charged K132 side chain is close to the images. The images 6 to 13 are stabilized by the negatively charged E136 and E140. The high barrier that appeared at around image 23 is attributed to steric hindrance by N82 and Coulomb repulsion by K154. In the case of route B, the conspicuous barrier is observed around images 23 and 24, which is produced by the steric hindrance by K154 and P273. Therefore, route B may be an unfavorable path for Na + conduction of this ChR structure.  case. This is because the barrier is located between two images in this case. The computational cost for � = 49 is larger than for � = 24. However, because the computational cost for the string method in the RismPath/SM is minor, it is recommended to employ a large number of image points.

Summary
In this letter, we have proposed a new method for finding the MFEP of ions and small molecule transportation through a channel protein based on the 3D-RISM theory combined with the string method, which is referred to as the RismPath/SM. Thanks to the 3D-RISM theory for evaluating the PMF of ions, the difficulties with sampling the ion and solvent distributions can be overcome. In the present study, we employed the zero-temperature string method for simplicity. It is noted that the finite temperature string method is also available, and it may be more effective in the case of a complicated PMF surface.
We applied the RismPath/SM to evaluate the MFEP of the Na + conduction pathway of the ChR closed state as an example. The method presented here is focused on the improvement of solvent distribution sampling for a given fixed protein structure. As expected, the structural fluctuation of the channel protein plays an important role in molecular transportation. Therefore, a combination with a method that can take the structural fluctuation of a protein into account, such as MD simulation, is strongly desired for highly accurate analysis. The self-consistent type of hybrid method of the 3D-RISM theory and MD simulation is suitable for this purpose [32,33].
In the method, the protein (and membrane, if any) structure is treated by MD simulation whereas the distributions of the ions and solvent molecules are evaluated using 3D-RISM theory. The combination of the MD/3D-RISM method with the efficient sampling methods of the protein structure is useful for providing the PMF surface to the RismPath/SM. Such studies with the RismPath/SM are in progress.