Globins are small respiratory proteins (140-160 residues) containing heme prosthetic group. This group is Fe-protoporphyrin IX ring. Heme is hosted in protein pocket defined by distal E- and proximal F-helices. Gaseous ligand (e.g. dioxygen) is bound reversibly by heme-Fe ion on distal side. On proximal side Fe atom is coordinated by a characteristic HisF8 residue. Despite large differences in their primary and quaternary structures, globins present a well-conserved tertiary structure organized in a three-on-three a-helical sandwich, so-called globin fold [6,7] (Fig. 1). Most globins serve either as O2 and CO2 transport proteins (hemoglobins) or as oxygen storage in muscles (myoglobins). Recent research extended the list of vertebrate globins by two new member: Neuroglobin (Ngb) and cytoglobin (Cgb) [8,9]. Both belong to the class of hexa-coordinated globins. Ngb is expressed in neural tissue and Cgb is expressed in all tissues at low concentration. Here is described ligand diffusion in nerve tissue minihemoglobin from Milky ribbon-worm (Cerebratulus lacteus) (CerHb) - another member of the neural globin family. Its sequence contains only 109 amino acids. This deletion of the A-helix and a shortening the H-helix in comparison with "standard" myoglobins may be a structural basis of the different affinity of CerHb to gaseous ligands than that measured for the other known globins. The structure of CerHb was discovered by Pesce A. at al. [1]. These authors indicated a wide channel leading from the heme cavity to the solvent but there is no useful information about a detailed route of the ligand from the cytosole to the heme pocked. No molecular dynamics simulations of ligand diffusion in this interesting system have been performed yet. Since it is important to know what residues affect ligand diffusion in this model we have initiated extensive, nanosecond scale MD studies of O2 diffusion in the CerHb matrix using CHARMM27 force field. It would be interesting to know what residues affect the ligand dynamics. Our results may serve as a good hint for making future kinetic experiments on CerHb variants

Structure of Cerebratulus lacteus minihemoglobin was obtained from Protein Data Bank (PDB). The id of this protein is 1KR7 (oxy form). The protein has 1676 atoms. Simulations of the minihemoglobin using the locally enhanced sampling (LES) [2,3] algorithm were carried out using the NAMD program [4] with CHARMM force field. Four, 1 nanoseconds simulations each with 1, 5, 10, and 15 copies of the dioxygen ligand were obtained. All simulations were performed in a model water box (TIP3, 49.5x49.0x50.0 A). The cutoffs for van der Waals and electrostatic interactions for all systems studied were 12 A. Switch distance was 8 A, and 1-4 scaling was used. The structures were initially frozen and equilibration of water molecules for 50 ps was performed. The heating was done during the next 100 ps of dynamics. In main simulations temperature was held on 300 K by Langevin dynamics (Langevin damping factor was 5). Final structures used in simulations had over 11 000 atoms. The analysis was performed using the VMD code [5]. Each trajectory was calculated on a cluster of 4xPentiumIV 2.0 GHz, 512 MB RAM running Linux Fedora 2.0 operation system. The Figures 7-11 were prepared using trajectory path script from VMD web site (authors: Andrew Dalke and Axel Kohlmeyer). Aminoacid composition of pathways for ligands diffusion were obtained using a script written by S.O. It was assumed that protein-ligand contacts (collisions) occurred when a distance between the ligand and an atom from hemogobin was lower than 2.5 A. The idea of the LES method is rather simple [2]. The whole system is divided on two subsystems: a big one and a small one. For example, in our case: The protein (CerHb - 1674 atoms)+water box (over 9000 atoms) and the ligand (dioxygen - 2 atoms). The several non-interacting copies of small subsystem (copies of the ligand) are build. These atoms feel a nonbonded (van der Waals and electrostatic) potential from the protein. On the other hand the big system (protein+water) experience an average potential from all these copies of the ligand. In this way enhanced atoms can occupy the same space, sampling of the conformational space is more efficient than in the case of just one small ligand and barriers are reduced [10]. This trick increase not only the sampling but also the transition rates for ligands. The LES trajectories with 10 copies may in principle give information equivalent to data obtained from an order of magnitude longer classical MD simulations. Some formal details of the LES method are given in Appendix.

In all trajectories Ca RMS (root-mean-square distance) with respect to the pdb protein structure do not exceed 1.5 A, thus we consider these data as reasonably stable (see Figs. 3-6). The lack of the A helix, a shorter H helix and other features indicate that a novel diffusion pathways in CerHb may be expected in comparison with a standard myoglobin picture [12,13]. 1 ns normal MD trajectory (1 copy) was too short to observe dioxygen diffusion out the protein matrix. The ligand sampled heme cavity and visited sometimes the other cavity, analogous to socalled (Xe..) in Mb [12,13] - see Fig. 7. The situation with 5 LES copies of the oxygen did not give any exit paths either (Fig. 8). Only LES trajectories with 10 and 15 copies had some successful 2 and 4 exit paths, respectively. Thust six different routes were found in total. Analysis of MD trajectories and O₂ collisions shows that at least two pockets accessible to O2 are present inside CerHb except the heme distal pocket (Fig. 6-7). The main heme pocket is composed of Tyr11, Leu14, Phe15, Phe25, Gln44, Thr48 and Phe10 residues. The first pocket is rather small and is composed of Tyr11, Phe27, Leu35, Ala40, Tyr41 (CavI). The second pocked (CavII) is composed of Ile52, Ile56, Asp61, Leu65, Leu86, Leu98, Ile102, Ala101. There are two alternative routes for trajectory with 10 copies of the dioxygen:
T1. The ligand left heme pocket through Ala40 and Lys43.
T2. From heme pocket ligand moved through the gate Tyr11. Then left CerHb through a tunnel composed of Tyr41, Asp8, Ala45, Cys42, Val7, Ala4, Val49 (between B and E helices). The trajectory with 15 copies gave the following four possible paths for ligand diffusion:
T3. The ligand left the heme pocket through the barrier at Lys47, Tyr51, Leu65.
T4. The ligand left the heme pocket through the barrier at Tyr21, Ala17, Asn81, Hsd18.
T5. The ligand moved from the heme pocket to the second pocket (Cav II) passing by Phe10, Tyr48, Ile52, Leu86. From there the ligand has left protein close to Ile56, Leu98, Ala101.
T6. Like in the route no. 5(T5) the ligand moved from the heme pocket to the second pocket, then left through the neighborhood of Ala62, Ala101, Asp104.
In routes number 5 and 6 the ligand left the CerHb between the E/F loop and H helix. As mentioned before, in 1 ns simulations with 1 and 5 copies dioxygen did not leave the protein.

In CerHb a very long (>10 ns) simulation will be necessary to obtain a classical diffusion path for a single dioxygen ligand at ambient temperature. The LES method helps to find out possible diffusion paths. Our data (10 and 15 LES ligands) suggest a possibility of multiple diffusion routes in CerHb. At least 6 paths were identified, The main barriers on crossing from the heme pocked to the next stages of ligand transport are Tyr11 and Phe10. In our opinion MD simulations give better insight into quite dynamic phenomenon of ligand diffusion in biopolymer interior than even careful, but static, inspection of the rigid protein X-ray structure.

The LES method [2] is based on the assumption that phase-space density can be written as product: