A molecular dynamics study on the diffusion and imprint ability of spectinomycin under different sizes of aniline oligomers

Molecularly imprinted polymers (MIP) are the polymers created by molecular imprinting techniques that leave cavities for the specific interactions with a template molecule and have been applied in molecular selectivity tasks. In this study, the molecular dynamics (MD) simulation technique was used to demonstrate that aniline oligomer could be developed as a potential MIP for detection and separation of the spectinomycin drug molecule for gonorrhea treatment. MD simulations were performed for the systems of a spectinomycin within aniline oligomers of different sizes. The mean square displacement (MSD) and the diffusivity calculated from MD simulations showed that the diffusion coefficient was significantly dropped when the length of aniline oligomer was greater than two. The diffusion coefficient of spectinomycin became the lowest within aniline trimers, corresponded to the highest atomic distribution of MIP around the template. Then, the specific cavity in MIP systems with and without spectinomycin was calculated to assess the stability of the cavity created by the template. The volume of a cavity created within the trimer system was closest to the spectinomycin volume and therefore became the optimal oligomer size for further development of MIP.


Introduction
Molecularly imprinted polymers (MIP) are polymeric materials created for a precise molecular recognition task by molecular imprinting techniques that leaves cavity for the specific interactions with a template molecule. They have been developed for many purposes, e.g., treatments, chromatography, and biosensing [1][2][3][4][5][6][7]. For pharmaceutical applications, selectivity of MIP for drug molecules has been an important issue [8][9][10]. Meanwhile, different types of MIP have been used to separate and remove the selected toxic compounds from the environment [11]. In order to design an MIP system to recognize a template molecule, matching up a wide variation of monomers and cross-linker molecules with the template of interest has become a challenge. Synthesizing all the subcomponent variants is extremely timeconsuming and required a large budget. Therefore, molecular simulation techniques and protocols were developed to overcome these challenges. Ab initio approaches, e.g., density functional theory, were used to investigate the interactions between monomer units of MIP and the template molecule of interest [12][13][14]. To consider more about thermodynamic effects and transport properties of polymers and templates, parameters obtained from DFT calculations [15] were transferred to the atomistic molecular dynamic (MD) simulations. In a study by Luo et al. [16] on the simulated MIP system made from methacrylate to detect a cholesterol molecule, the ability of MIP to immobilize the template was measured through the mean-square displacement (MSD) and diffusion coefficients of the constituent molecules. Both methods provide the detail of molecular interactions, which are essential for the further designs of MIP.
In this study, we tried to develop an atomistic MD simulation and a trajectory analysis protocol for MIP systems. A case study of the polyaniline(PANI)-based MIP built on the spectinomycin template was picked. Spectinomycin is an aminocyclitol antibiotic produced by Streptomyces spectabilis. The common use of spectinomycin is for the treatment of gonorrhoea in human. Veterinarians used the intramuscular injection of spectinomycin for cattles got infected from consuming contaminated feed. However, the spectinomycin residue in dairy products may cause an allergic skin reaction for humans [17]. Designing an MIP system suitable for the spectinomycin template could be developed into a highly selective sensing platform. Polyaniline (PANI) is a semiflexible conducting polymer previously used for such templates as chloramphenicol, ascorbic acid, and tryptophan [18][19][20]. Aniline monomers or oligomers become interacted with the template, and an MIP was formed after the polymerization process. From our simulations, the effects of size or degree of polymerization on the properties of aniline-based MIP were investigated. Mobility and diffusion of oligomers and the spectinomycin template were measured. Also, an analysis on the cavity created by the template molecule was performed to address the imprinting ability of the MIP model.

Starting structures and topology files of MIP models
Each of the four MIP models consisted of two components: (1) the spectinomycin template molecule and (2) the functional oligomers of degree of polymerizations 1 to 4. The structure of spectinomycin ( Fig. 1a.) was obtained from the PubChem database [21], and the structures of aniline oligomers (Fig. 1b) were built by the Avogadro software [22]. Then, topology files of all molecules were generated from the Automated Topology Builder (ATB) webserver, in which a DFT calculation was performed for each structure to obtain the partial charge distribution q i on each atom. The information on the covalent bonds, e.g., constants for stretching ( k i ), bending ( k ), and torsional motion ( k ), was then assigned according to the GROMOS 54A7 force-field parameter set [23], in which potential energy functions for covalent bonding were contributed for each atom i and its neighbors and were given in terms of atomic coordinates as follows: For the non-covalent terms in (2), van der Waals and electrostatic interactions were defined for atom pairs ( i, j ) through the sum of van der Waals radius ij of the pairs, the van der Waals interaction strength ij , and the partial charges q i and q j as follows:

Molecular dynamics (MD) simulations
For each calculation step t in an MD simulation, each degree of freedom i was subjected to the force F i that came from the gradient or the first derivative of forcefield potential U = U covalent + U noncovalent with respect to position ⃗ r i , shown as follows: Then, the numerical solution of the classical equation of motion ⃗ F i = m i 2 ⃗ r i t 2 to obtain the trajectory was solved by the velocity-Verlet algorithm, providing the updated position and velocity at the next time step t + Δt as follows: To prove that aniline oligomers could form an MIP to capture a spectinomycin, the approach similar to a previous study of Zhi-Jian et al. [24] was used. For each of the four systems, a 4nm × 4nm × 4nm simulation box was created to surround a spectinomycin molecule. Then, the four simulation boxes were filled with aniline oligomers (see Fig. 1c). After an energy minimization, 100-ns MD production runs under an NPT ensemble were carried out at 300 K and 1 atm. The time step of 2 fs was used along with the PLINC algorithm for holonomic constraints of each covalent bond involving a hydrogen atom. After that, the spectinomycin was removed from each system, and the simulation was continued for 100 ns under the same condition to observe the imprinting ability of template on the MIP. All simulations were performed by the GROMACS 5.1.2 package [25].

Analysis
To quantitatively monitor the formation of MIP, mean square displacement (MSD) of the monomers and spectinomycin was calculated from the trajectories as follows: Diffusivity ( D ) was then calculated from the slope of the between the first 10% and the last 10% , in which MSD linearly increased with time and the diffusivity is defined by the following: In addition of measuring diffusivity of both templates and MIP molecules from MSD, diffusivity can also be calculated from the Stokes-Einstein relation as follows: when the dynamic viscosity was calculated for the representative nitrogen atoms by using the viscosity calculation tool embedded within the "gmx energy" module. The viscosity was calculated from the first 10 ns of MD trajectories to avoid the turbulence flow. As our polymer systems were assumed to comprise of spherical monomer beads, the Stokes-Einstein relation was valid, and the radius r could be found from half the median distance between aniline monomers from different aniline oligomer chains.
Then, to quantitatively investigate the intermolecular interactions between the spectinomycin template and (4) surrounding aniline oligomers, e.g., the distribution of the nitrogen atoms within aniline oligomer as either the main hydrogen bond donors or acceptors about the hydrogen and oxygen atoms of spectinomycin, radial distribution functions (RDF; g(r)) defined by were calculated to visualize the pair correlation between these two groups of atoms separated by the distance r . The term dn r ∕4 r 2 dr represented the density of the n r acceptor atoms at the distance r from a donor atom. The density of the acceptor atoms within a very thin spherical shell with radius r centered at a donor atom was then normalized by the local density within the whole spherical volume of radius half the simulation box size. Further detail of hydrogen bonding and hydrophobic contacts of the final simulated snapshot was then visualized through the LIGPLOT program [26].
Finally, to assess the imprinting ability of MIP models, the cavity volume was measured through the flood-fill method [27]. Snapshots were taken from each simulation, and water molecules were added into the cavity created by spectinomycin template by the "gmx insert-molecules" command. Then, the cavity volume estimated through the maximum number of water molecules and the density of water at 300 K was compared for simulations with and without the spectinomycin template at different oligomer lengths. Figure 2 displayed the mean square displacement of spectinomycin (red line) within aniline oligomers along with the corresponding superimposed snapshots. MSD of aniline monomers and spectinomycin within the monomers were found significantly higher than MSD of other simulations. Diffusion coefficients of spectinomycin and aniline monomers were found at 123.000 × 10 −9 cm 2 ∕s and 2708.000 × 10 −9 cm 2 ∕s , respectively (Fig. 2a). Conformational snapshots of the spectinomycin demonstrated its high mobility, suggesting that anilines were in the liquid phase. In Fig. 2b, aniline monomers were assumed to be polymerized into dimers, and the matrix became glassy. The glassy state of aniline dimer was indicated by an abrupt decrease in the diffusion coefficient from 2708.000 × 10 −9 cm 2 ∕s of the monomers to 2.893 × 10 −9 cm 2 ∕s of the dimers (by 936.05-fold). This also resulted in the diffusion coefficient of spectinomycin to drop from 123.000 × 10 −9 cm 2 ∕s of the monomers to 8.000 × 10 −9 cm 2 ∕s of the dimers (by 15.37-fold). In Fig. 2c, as the oligomers became trimers, the diffusion coefficient of aniline trimers 2.000 × 10 −9 cm 2 ∕s remained less than that aniline dimers (by 4.00-fold). In the same way, the diffusion coefficient of spectinomycin became as low as 0.889 × 10 −9 cm 2 ∕s . For the case of aniline tetramers in Fig. 2d, the diffusion coefficient of aniline tetramers was found at 0.466 × 10 −9 cm 2 ∕s , suggesting stronger interactions within the oligomer matrix due to the increased number of covalent bonds. The diffusion coefficient of spectinomycin within the tetramers became 0.055 × 10 −9 cm 2 ∕s , and the whole system was completely in a glassy state. In the same way, when considering diffusivity calculated from the Stokes-Einstein equation shown in Table 1, the spectinomycin within aniline monomers was with a significantly high mobility due to a low dynamic viscosity. As the size of MIP molecules became larger with two or more aniline monomer units, the diffusion coefficient became significantly decreased, which also suggested the transition into a glassy state.

Interaction analysis
Apart from many hydrophobic contacts between the aromatic rings of aniline oligomers and the nonpolar functional groups of spectinomycin, specific interactions between the spectinomycin and the oligomers were characterized by the hydrogen bonds. Figure 3a displayed the radial distribution function (RDF) measured between the group containing oxygen atoms of spectinomycin as either  the h-bond donors or acceptors and the group containing nitrogen and hydrogen atoms of aniline oligomers. Meanwhile, Fig. 3b displayed the RDF measured between oxygen atoms of spectinomycin and the polar hydrogen atom of aniline oligomers. The sharp first RDF peaks observed for both plots in Fig. 3 suggested that the interactions between donor and acceptor atom groups were favorable, and hydrogen bonds were likely to form along with the hydrophobic contacts, as confirmed by the interaction analysis of the last snapshots after 100 ns at all oligomer lengths (Fig. 4). The last MD simulation snapshot for the system containing aniline monomers (Fig. 4a) displayed an only hydrogen bond along with hydrophobic contacts between the spectinomycin template and five aniline monomers. Surprisingly, as the aniline became dimeric (Fig. 4b), the number of hydrophobic contacts increased to seven, but the number of hydrogen bond became one.
For the system of aniline trimers (Fig. 4c), the number of hydrogen bonds increased from one to four, and the number of hydrophobic contacts was similar to the dimer system before decreasing to five contacts for the tetrameric systems (Fig. 4d). The second and the third RDF peaks indicated weak attractions between the second and the third nearest neighbors and indicated the order emerged within the systems of oligomer chains that contained aromatic rings. For the case of aniline monomer, the second and third RDF peaks almost disappeared, and the RDF value was converged to 1 as the distance between atoms became larger than 1 nm, signifying the disordered liquid phase of aniline monomers under the simulated conditions. The RDF between h-bond donors and acceptors for all the polymerized anilines in dimeric (red line), trimeric (green line), and tetrameric (black line) displayed the glassy behavior due to the formation of oligomer networks, interacting through pi stacking. Small correlations were found at further distance than 1 nm, indicating that stronger interactions between oligomers caused an amorphous matrix. The highest RDF peak for trimers and tetramer suggested the highest propensity for hydrogen bonding and was in concurrence with the lowest mobility of spectinomycin.

Cavity volume
When spectinomycin was removed, a cavity appeared within the matrix and should disappeared shortly within an NPT simulation if the matrix was in the liquid state with fast relaxation time. However, for an MIP with good imprinting ability, volume of the created cavity should be maintained. From our MD simulations of MIP systems with four different oligomer sizes, cavity volume calculated by the flood-fill method was used to measure the cavity maintenance after template removal.
As the system with aniline monomers was in liquid state, the fast relaxation time of aniline promptly closed the cavity after the template removal and indicated that aniline monomer was unable to provide the matrix for an MIP. On the other hand, the cavities within the aniline dimers, trimers, and tetramers created by the spectinomycin template were maintained, as shown in Fig. 5, where a number of water molecules was inserted into the cavity by the flood-fill method both before and after the spectinomycin removal, followed by a 50-ns equilibration.
After that, the cavity volumes within the aniline dimer (Fig. 5a), trimer (Fig. 5b), and tetramer (Fig. 5c) systems were calculated from the maximum number of water molecules presented in the cavity. The calculated volume for each oligomer systems after 100 ns without the template was compared to when Fig. 3 Radial distribution functions of the a group containing hydrogen and nitrogen atoms of aniline monomer (blue), dimer (red), trimer (green), and tetramer (black) relative to oxygen atoms of spectinomycin and b the polar hydrogen atoms within aniline monomer (blue), dimer (red), trimer (green), and tetramer (black) relative to oxygen atoms of spectinomycin the template was presented, as the swelling of MIP was shown to be important for template detection [28]. Table 2 showed that the cavity volume of small aniline oligomer of lengths 2-4 monomers was in the range of 1.5-2.5 times the volume ~ 600 cm 3 of spectinomycin. Cavity volume of the aniline dimer MIP model was shown to increase by twofold. This volume expansion observed for the dimeric system provided a clue on the strength of cohesive forces between the bulk aniline dimers. The trimeric and tetrameric systems were with their cavity volumes decreased after the spectinomycin removal. Compared to other systems, aniline trimers formed a cavity with the volume closest to that of spectinomycin, in concurrence with the highest degree of hydrogen bonding and hydrophobic contacts with the spectinomycin template and the highest first RDF peaks. Moreover, the smallest volume difference of the trimeric system after the spectinomycin removal indicated the greatest imprinting ability. From all these results, aniline trimer possessed the optimal degree of polymerization for creating an MIP for spectinomycin detection.

Conclusions
In this study, an atomistic MD simulation and a trajectory analysis protocol for MIP systems were developed. The polyaniline(PANI)-based MIP built on the spectinomycin Fig. 4 The interaction of spectinomycin with a aniline monomer, b aniline dimer, c aniline trimer, and a aniline tetramer determined by LIG-PLOT [26]. Hydrogen bonds are represented by green dashed lines, and hydrophobic interactions are represented by red arcs template was picked as a case study. Our results showed that the polymerized aniline, even as dimers, displayed a glassy behavior. Moreover, it was found that aniline oligomers specifically interacted with the spectinomycin template through hydrogen bonds, and the aniline trimer could maintain the cavity after the template removal and was fit to the volume of spectinomycin template, proving as a good candidate for developing the MIP for spectinomycin and sensing platforms for food industries to prevent side effects from the consumption of spectinomycin residues.

Data availability
The authors confirm that the data supporting the findings of this study are available within the article and its supplementary material.
Code availability Not applicable.

Declarations
Conflict of interest All authors have participated in (a) conception and design or analysis and interpretation of the data, (b) drafting the article or revising it critically for important intellectual content, and (c) approval of the final version. This manuscript has not been submitted to, Fig. 5 Cavity created by the spectinomycin molecule within the MIP models of a aniline dimer, b trimer, and c tetramer before and after the spectinomycin removal. Aggregated aniline oligomers were represented in blue, while water molecules filled within the cavities were represented in gray