.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "examples/gallery/structure/modeling/water_exchange_noexec.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_examples_gallery_structure_modeling_water_exchange_noexec.py>`
        to download the full example code

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_examples_gallery_structure_modeling_water_exchange_noexec.py:


Cavity solvation in different states of an ion channel
======================================================

In selective sodium/potassium channels, the internal cavity of the pore
is walled off from the solvent if the channel is closed.
Upon activation, the internal gate opens and exchange of water molecules
between the cavity and the bulk medium is possible.

Therefore, one can track the exchange rate of water molecules between
the cavity and bulk to evaluate if a pore is open, closed, or in a
transition between the two. Here, we used the distance between water
molecules and residues located in the central cavity to evaluate if
persistant water exchange takes place in different structures of the
HCN4 channel.

The trajectories and template structure are not included in this
example.
However, the trajectories are based of publicly accessible structures
of the open (PDB: 7NP3) and closed (PDB: 7NP4) state.

.. image:: ../../../scripts/structure/modeling/water_exchange.png

.. GENERATED FROM PYTHON SOURCE LINES 24-122

.. code-block:: Python


    # Code source: Daniel Bauer, Patrick Kunzmann
    # License: BSD 3 clause

    import matplotlib.pyplot as plt
    import numpy as np
    from pylab import polyfit
    import biotite
    import biotite.structure as struct
    import biotite.structure.io.gro as gro
    import biotite.structure.io.xtc as xtc


    def water_in_prox(atoms, sele, cutoff):
        """
        Get the atom indices of water oxygen atoms that are in vicinity of
        the selected atoms.
        """
        cell_list = struct.CellList(atoms, cell_size=5, selection=atoms.atom_name == "OW")
        adjacent_atoms = cell_list.get_atoms(atoms[sele].coord, cutoff)
        adjacent_atoms = np.unique(adjacent_atoms.flatten())
        adjacent_atoms = adjacent_atoms[adjacent_atoms > 0]
        return adjacent_atoms


    def cum_water_in_pore(traj, cutoff=6, key_residues=(507, 511)):
        """
        Calculate the cumulative number of water molecules visiting the
        pore.
        """
        protein_sele = np.isin(traj.res_id, key_residues) & ~np.isin(
            traj.atom_name, ["N", "O", "CA", "C"]
        )
        water_count = np.zeros(traj.shape[0])
        prev_counted_indices = []
        for idx, frame in enumerate(traj):
            indices = water_in_prox(frame, protein_sele, cutoff)
            count = (~np.isin(indices, prev_counted_indices)).sum()
            if idx != 0:
                count += water_count[idx - 1]
            water_count[idx] = count
            prev_counted_indices = indices
        return water_count


    # Calculate the cumulative number water molecules visiting the pore
    # for the open and closed state
    counts = []
    for name in ["apo", "holo"]:
        gro_file = gro.GROFile.read(f"{name}.gro")
        template = gro_file.get_structure(model=1)
        # Represent the water molecules by the oxygen atom
        filter_indices = np.where(
            struct.filter_amino_acids(template) | (template.atom_name == "OW")
        )[0]
        xtc_file = xtc.XTCFile.read(f"{name}.xtc", atom_i=filter_indices)
        traj = xtc_file.get_structure(template[filter_indices])
        cum_count = cum_water_in_pore(traj)
        counts.append(cum_count)
    time = np.arange(len(counts[0])) * 40 / 1000


    # Linear fitting
    open_fit = polyfit(time, counts[0], 1)
    closed_fit = polyfit(time, counts[1], 1)

    fig, ax = plt.subplots(figsize=(8.0, 4.0))
    ax.plot(time, counts[0], label="open pore", color=biotite.colors["dimgreen"])
    ax.plot(
        time, open_fit[0] * time + open_fit[1], linestyle="--", color="black", zorder=-1
    )
    ax.plot(time, counts[1], label="closed pore", color=biotite.colors["lightorange"])
    ax.plot(
        time, closed_fit[0] * time + closed_fit[1], linestyle="--", color="black", zorder=-1
    )
    ax.set(
        xlabel="Time / ns",
        ylabel="Count",
        title="Cumulative count\nof individual water molecules visiting the pore",
    )
    ax.legend()
    ax.annotate(
        f"{open_fit[0]:.1f} per ns",
        xy=(20, 20 * open_fit[0] + open_fit[1] + 100),
        xytext=(20 - 5, 20 * open_fit[0] + open_fit[1] + 1300),
        arrowprops=dict(facecolor=biotite.colors["darkgreen"]),
        va="center",
    )
    ax.annotate(
        f"{closed_fit[0]:.1f} per ns",
        xy=(30, 20 * closed_fit[0] + closed_fit[1] + 100),
        xytext=(30 + 2, 20 * closed_fit[0] + closed_fit[1] + 1300),
        arrowprops=dict(facecolor=biotite.colors["orange"]),
        va="center",
    )
    fig.savefig("water_exchange.png", bbox_inches="tight")

    plt.show()


.. _sphx_glr_download_examples_gallery_structure_modeling_water_exchange_noexec.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: water_exchange_noexec.ipynb <water_exchange_noexec.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: water_exchange_noexec.py <water_exchange_noexec.py>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_