Simulation of benzylpenicillin molecule distribution in slit-shaped Si nanopores

A molecular dynamics study of the behavior of benzylpenicillin molecules in slit-shaped nanopores was carried out. A model silicon material with a pore size from 10 to 50 nm was chosen as a nanoporous structure. The interaction between benzylpenicillin molecules was described by a pair potential, built on the basis of modelling the molecule behavior by all-atom force fields. It was shown that an adsorbed layer of benzylpenicillin molecules is formed near the pore walls. With a decrease in the pore size, the maximum density of molecules in the adsorbed layer decreases, while the fraction of adsorbed molecules in the whole pore increases.


Introduction
The behavior of gas and liquid in nanoscale pores may drastically differ from their evolution in a free state [1]. It was established experimentally that a decrease in the size of nanopore leads to a shift of the coexistence curves and a decrease of critical points [2]. In particular, the 3D behavior of the system changes to 2D as the size of carbon and hydroxyapatite slit-shaped pores decreases [3,4]. Reduction of the diameter of cylindrical nanopores in some zeolites and SiO2 leads to 1D gas behavior in them [5,6]. The change in the behavior of gas and liquid in nanoscale pores compared with their behavior in the free state is largely due to the effect of adsorption [7]. The effect of increase in the volume fraction of gas or liquid molecules in the near-wall region at decreasing nanopore size enhances the influence of the adsorption on the properties and behavior of the substance inside the nanopore [8,9]. In turn, a fluid in pores of different sizes affects the mechanical properties of porous materials in different ways [10].
Studying the characteristics of the adsorption processes in porous materials would deepen the existing knowledge about the features of substance behavior inside them, which is of considerable interest for the creation of scientific foundations for the development of catalysts, devices for gas sensing and nanoscale drug carriers for targeted delivery to affected tissues. In the development of various practical applications, scientific knowledge is needed about the phenomena associated with surface effects due to the interaction between the substance and pore walls, the size, morphology and distribution of pores [11].
In connection with the foregoing, in this work, a computer simulation was performed to study the effect of the silicon nanopore size on the features of benzylpenicillin molecules adsorption.

Methods
The peculiarities of the drug molecule adsorption in the slit-shaped nanopores were studied within the framework of the molecular dynamics approach [12][13][14][15][16]. The distance between the nanopore walls ranged from 10 to 50 nm. In present simulations the benzylpenicillin molecule (PCN) treated as a material point that interacts with the silicon surface and with other PCN molecules according the tabulated potential, built on the basis of calculating the change in the free energy of PCN-Si and PCN-PCN interactions as functions of distance. To calculate the energy of interaction of benzylpenicillin with a regular silicon surface at the molecular level, using all-atom force field, a model of benzylpenicillin in a state with a deprotonated carboxyl group near the (100) Si surface with fully hydrated water (TIP3P model [17]) with a counterion of sodium was prepared. To estimate the potential of the interaction of benzylpenicillin molecules with the silicon surface, the computational cell contained 4473 atoms, including 1216 water molecules. The simulation time was 260 ns with an integration step of the equations of motion 1 fs. In the calculations for estimating the pair PCN-PCN potential, the volume of the hydrated computational region was V=57.262 nm 3 , the model contained two benzylpenicillin molecules, two Na + ions and 1939 water molecules (a total of 5899 atoms). The free-energy change profile of PCN-Si and PCN-PCN interaction as a function of distance was estimated from a series of bidirectional steered molecular dynamics simulations with a mean-force potential analysis.
All simulations were carried out in the LAMMPS package [18] at a constant temperature of 310 K. Simulation of the behavior of drug molecules in silicon nanopores was carried out for densities from 0.035 to 0.35 molecules/nm 3 . The equilibrium distribution of PCN molecules over the cross section of the simulated pore was calculated on the basis of averaging the simulation results over a time interval of 100 ns. Visualization of the simulation results was carried out in the OVITO package [19].

Results and discussion
The simulations showed that the nanopore size has a significant effect on both the time needed to reach an equilibrium state and the features of the distribution of PCN molecules in the nanopore. In the initial state, PCN molecules were uniformly distributed over the nanopore volume (Fig. 1a). In the process of the simulated system relaxation, the molecules redistribute over the cross section of the nanopore (Fig. 1b, c).
With the same density of molecules in larger nanopores, the relaxation time of the simulated system increases approximately linearly (Fig. 2a). In the process of reaching the equilibrium distribution of molecules, the adsorbed layer is formed (Fig. 1c), which leads to a change in the interaction energy between the wall of the nanopore and PCN molecules (Fig. 2b). The density of molecules in the adsorbed layer is controlled by the peculiarities of: 1) PCN intermolecular interaction, and 2) interaction of PCN molecules with the wall of nanopore. The balance of these interactions determines the molecular density in the adsorbed layer. Note that the higher adsorbed molecules density corresponds to the larger average repulsive force between the molecules in the layer. In the case of PCN density of 0.35 molecules/nm 3 , the maximum attraction of the PCN molecule to the wall is at a distance of 0.45 nm, which corresponds to the coordinate of the maxima on the density distribution curves over the nanopore cross section. The density distribution of PCN molecules near the nanopore wall after reaching thermodynamic equilibrium is shown in Fig. 3a. It can be seen that the density of molecules in the adsorbed layer increases with increasing width of the slit-shaped nanopore. The maximum local density of molecules in the adsorbed layer corresponds to a nanopore with 50 nm width and is approximately 7.5 molecules/nm 3 , and the minimum density is 5.3 molecules/nm 3 for the nanopore with 10 nm width. A decrease in the average density of molecules in nanopores does not affect the features of their distribution over the cross section for different nanopore widths.
Simulations showed that an increase in the width of the nanopore and the average total density of PCN molecules contained in it leads to a significant increase in the maximum density of molecules in the adsorbed layer. At the same time, the fraction of molecules in the adsorbed layer decreases with increasing width of the nanopore for all studied values of the average density of PCN (Fig. 3b).

Conclusions
Our molecular dynamics simulations have shown that the evolution of PCN molecules in porous silicon is substantially different compared to their evolution in free state. An adsorbed layer of PCN molecules is formed in the region of the pore walls, the density of which may be several times higher than the average density of molecules in the bulk. It was found that an increase in the width of the slit-shaped nanopores leads to an increase in the density of molecules in the adsorbed layer. With a decrease in the size of the slit-shaped nanopore, the fraction of adsorbed molecules increases.
The work was financially supported by the Russian Federation via the Ministry of Science and Higher Education of the Russian Federation (project identifier RFMEFI60718X0202).