Molecular dynamics performance for coronavirus simulation by C, N, O, and S atoms implementation dreiding force field: drug delivery atomic interaction in contact with metallic Fe, Al, and steel.

Coronavirus causes some illnesses to include cold, COVID-19, MERS, and SARS. This virus can be transmitted through contact with different atomic matrix between humans. So, this atomic is essential in medical cases. In this work, we describe the atomic manner of this virus in contact with various metallic matrix such as Fe, Al, and steel with equilibrium molecular dynamic method. For this purpose, we reported physical properties such as temperature, total energy, distance and angle of structures, mutual energy, and volume variation of coronavirus. In this approach, coronavirus is precisely simulated by O, C, S, and N atoms and they are implemented dreiding force field. Our simulation shows that virus interaction with steel matrix causes the maximum removing of the virus from the surfaces. After 1 ns, the atomic distance between these two structures increases from 45 to 75 Å. Furthermore, the volume of coronavirus 14.62% increases after interaction with steel matrix. This atomic manner shows that coronavirus removes and destroyed with steel surface, and this metallic structure can be a promising material for use in medical applications. © OWZ 2020.

Introduction

Coronaviruses are a group of connected RNA viruses that origin illnesses in birds and mammals. In humans, these viruses cause some problem in respiratory system infections that can range from mild to lethal. Mild diseases contain some cases of cold and fatal varieties can cause diseases such as MERS, COVID-19, and SARS. There are any vaccines or antiviral drugs to stop or remedy human coronaviruses infections yet. Coronaviruses constitute the subfamily Orthocoronavirinae, in the family Coronaviridae, order Nidovirales, and realm Riboviria [1, 2]. Historically, coronavirus was recognized in the 1930s for the first time [3]. This virus depicted in Fig. 1.

Fig. 1

a The common structure of coronavirus recognized by doing some experiment researches. b All-atom depiction of the coronavirus simulated using the MD approach [4]

Researchers reported a new respiratory infection of chickens in North Dakota in 1931. The infection of chicks was distinguished by listlessness and gasping. The death rate of the chickens was up to 90% [5, 6]. Transmissible gastroenteritis and mouse hepatitis virus were detected in the 1940s [7]. Chemically, these structures are enveloped viruses with a nucleocapsid of helical symmetry and a single-stranded RNA genome which located in an icosahedral protein shell [8]. This atomic structure is one of the biggest thorugh RNA viruses. The genome size of coronaviruses ranges from about 26 to 32 kilobases [4, 9]. Today, numerous researches have been done on the characteristics of this virus. Masters addressed the present state of comprehension of coronavirus genome packaging, which has mainly been studied in 2 prototype species, transmissible gastroenteritis virus, and mouse hepatitis virus [10]. Lu er al. [11] phylogenetic analysis shows that bats, an animal sold at the seafood market in Wuhan, might be the original host of this virus in human.

Furthermore, the structural analysis of this research group suggested that coronavirus might be able to bind to the angiotensin-converting enzyme 2 receptors in humans. Then researchers introduce evolutionary features of the coronaviruses and successfully predicted new coronavirus outbreak by pointing out that novel pathogenic variants will readily appear from very various severe acute respiratory syndrome-related coronaviruses of the bat origin through their close coexistence and high genetic recombination potential [12]. Contrary to many studies, it is not known how the coronavirus interacts with the human environment, such as metallic structures. In this work, we study the coronavirus atomic interaction with metallic surfaces with MD [13-18]. Investigations are provided for the coronavirus and metallic matrix (Fe, Al, and steel structures) interactions at various temperatures (300–350 K) in this paper. For this purpose temperature, potential energy, center of mass (COM) distance, COM angles, mutual energy, and volume of simulated structures are reported.

Numerical method

In this computational study, we used MD simulation to calculate the atomic behavior of coronavirus in contact with metallic surfaces such as Fe, Al, and steel. Molecular dynamics method is a way for calculating the dynamical evolution of atoms and molecules. In MD method, atoms and molecules are allowed to atomic interaction and give a view of the position and velocity changes of the total system. In the research, all MD simulations were done via LAMMPS [19-22]. This computational package designing began in the 1990s by Sandia and LLNL laboratories. To use this package to calculate the atomic behavior of coronavirus, this atomic structure, H2O molecules, and metallic matrix were simulated. For this, coronavirus structure was situated in the middle region of MD simulation package, while the other regions of MD simulation box was filled by H2O structure. This initial atomic positions prepared via Packmol package [23]. Figure 2 displays the MD simulation package at the front, top, and perspective views, which were depicted by open visualization tool (OVITO) [24]. Particle base methods such as MD, LBM [25-36], or nanoparticles dispersed in the base fluid can be estimated as the various approaches of investigation through such that problem [37-51].

Fig. 2

Schematic of coronavirus/H2O molecules mixture and metallic matrix (Fe structure) which prepared via Packmol box at a front, b top, and c perspective views

In this molecular dynamics simulations, periodic boundary conditions were used in x, y, and z directions. For preparing the initial temperature of simulated structures, the thermostat was used in the MD simulation box to equilibrate this physical parameter at 300 K, 325 K and 350 K with 1 femtosecond time step. After 1,000,000 MD simulation time steps, at the selected temperature, the simulated atomic structures were equilibrated; then, micro-canonical ensemble (NVE) was implemented to describe the virus and metallic surfaces atomic interaction. To simulate the coronavirus atomic structure, we use the dreiding force field [52]. This interatomic force field is the appropriate choice to biological mixtures MD simulation. In dreiding force field, the potential of various atoms was described as a non-bonded and bonded forces. Non-bond force between particles in dreiding force field described by the Lennard–Jones (LJ) formalism. This formalism is a mathematically simple relation about the interatomic force between a pair of particles. This simple relation was stated by Lennard Jones [53] for the first time. This simple atomic interaction stated as Eq. (1):

In Eq. (1), σ is the distance at which the function is zero; ε is the depth of the potential well, and r is the distance between the two atoms. In MD simulations, both σ and ε constants related to the kind of the particles in the MD simulation package. The length scale parameter, cutoff radius, and energy for various atoms in coronavirus simulation are classified in Table 1 [52].

Table 1

The ε and σ constants for LJ interaction in coronavirus MD simulation [52]

	ε (kcal/mol)	σ (A)
Carbon	0.3050	4.180
Nitrogen	0.4150	3.995
Oxygen	0.4150	3.710
Sulfur	0.3050	4.240

The bonded forces consist of bond angle bend, bond strength, and dihedral angle torsion terms. The bond and angle strength stretch in dreiding force field are calculated by harmonic oscillator equations as (2) and (3) formulas:

In these formulas, K and K are harmonic oscillator constants. θ0 is the equilibrium value of angles, and r0 is the atomic bond length. K and K constants in coronavirus MD study were selected at 100 (kcal/mol)/degree2 and 300 (kcal/mol)/A2, respectively. Furthermore, the r0 and θ0 constants in our MD simulations are classified in Table 2 [52]. Dihedral term in atomic interaction described with Eq. (4) and its coefficients are chosen from a dreiding force field [52]:

Table 2

The r0 and θ0 rates of coronavirus structure in our MD simulations [52]

Bond/angle	r0 (Å)	θ0 (°)
CC	1.530	–
NC	1.462	–
OC	1.420	–
SC	1.800	–
ON	1.352	–
CCC	–	109.471
CNC	–	106.700
CSC	–	92.100
NCC	–	109.471
NCN	–	109.471
OCC	–	109.471
OCO	–	109.471
SCC	–	109.471
NOC	–	104.510
ONC	–	106.700

In Eq. (4), K is an oscillator constant with + 1 or − 1 rates and the integer number is n [52].

Today, various atomic models such as SPC, TIP4P, and TIP3P are used to molecular dynamics simulation of water atomic structure. In SPC model, three sites were used for the electrostatic interactions and the positive charges on the hydrogen atoms were neutralized by a negative charge. The interaction between H2O molecules were simulated using a LJ potential. These interaction parameters for atomic water structure are reported in [54, 55]:

O weight = 15.9994 u

H weight = 1.008 u

O atomic charge = − 0.820

H atomic charge = 0.410

ε constant for O–O interaction = 0.1553 kcal/mol

σ constant for O–O interaction = 3.166 Å

ε and σ constants for O–H and H–H interactions = 0.0

r0 constant for OH bond = 1.0 Å

θ0 constant for HOH angle = 109.47◦

Furthermore, metallic structures in our MD simulations described by inserted atom model (EAM) force field [56]. EAM force field is represented by Eq. (5):

In this equation, r is the distance between i and j particles, ϕ is a pair-wise potential function, ρ is the contribution to the electron density from particle j of type β at the location of particle i, and F is an embedding function that stats the energy required to place particle i of type α into the electron distribution. After atomic modeling and force field implementing, the atomic manner of coronavirus in the vicinity of metallic matrix estimated. For computations of these atomic evolutions, Newton’s second law [13, 14],

Which leads to [13, 14]:

Temperature from the Gaussian distribution [13, 14]:

Association of Eq. (6) is provided by velocity Verlet algorithm.

In these Eqs. (9) and (10), , is the final coordination and velocity of particles and v(t) and r(t) are the rate of these physical parameters at t = 0. Finally, we can say that MD simulations of this work in 2 steps [57-59]:

Step A Coronavirus and metallic matrix were simulated at 300, 325, and 350 K with 1 fs time step.

Step B In the second step of our molecular dynamics simulation, the atomic interaction between the virus and metallic matrix was carried out for 1,000,000 MD simulation time steps.

Results

Equilibration process of atomic structures

In the first step of this molecular dynamics simulations, the atomic structure of the coronavirus and metal surfaces was studied. Figure 3 shows the atomic arrangement of structures in our simulation box [52]. Numerically, structures atomic stability are described by reporting of temperature and potential energy of them at 300 K, 325 K, and 350 K. Temperature variation of various structures is depicted in Fig. 3. From this figure, we can say that all structures thermally equilibrated after 1,000,000 MD simulation time steps passing. Figure 4 displays the potential energy. We can say that the simulated structures potential energy converged after 1,000,000-time steps. This result shows that dreiding force field has good ability in bio structures simulations. Numerically, the potential energy of the coronavirus structure with Fe matrix decreases from − 28,545 to − 25,665 eV by temperature increasing from 300 to 350 K. This atomic manner shows the kinetic energy increasing in simulated structures which this parameter increasing cause the mean distance of atomic structure rises. It can be said that simulated structures with Fe matrix have better stability rather than Al and Steel matrix and by temperature increasing this stability decreases (see Fig. 4). Furthermore, coronavirus atomic stability reaches to minimum rate in the presence of steel matrix, while this atomic behavior of the coronavirus makes the steel matrix work best in medical applications.

Fig. 3

Temperature variation of simulated structures at a various initial temperatures and b various metallic matrix as a function of MD simulation time steps

Fig. 4

Potential energy variation of simulated structures at a various initial temperatures and b various metallic matrix

Dynamical evolution

The COM of the atomic structures is one of the geometrical. The interatomic force between the coronavirus and metallic surface is repulsive (Fig. 5). Numerically, the distance of coronavirus and Fe atoms varies from 45 to 57 Å at 300 K. This physical parameter increases by temperature increasing from 300 to 350 K (see Fig. 6 and Table 3). As reported in Table 3, steel and aluminum matrix interaction with coronavirus shows a similar manner and steel matrix shows maximum atomic distance from the coronavirus after 1,000,000-time steps. Furthermore, as seen in Table 4, the angle of the coronavirus with all matrix surface fluctuates with temperature increasing. This geometrical parameter has 92°, 95°, and 88° at 300 K, 325 K and 350 K for Fe matrix surface, respectively. So by increasing temperature, the angle of the coronavirus does not follow a logical relation.

Fig. 5

Time evolution of coronavirus and metallic matrix with COM distance for a 0.5 ns, b 1 ns, c 1.5 ns, and d 2 ns

Fig. 6

COM distance variation of coronavirus and Fe matrix with MD simulation time steps and sample temperature

Table 3

COM distance of coronavirus and metallic matrix

Temperature (K)	Fe atomic distance (Å)	Al atomic distance (Å)	Steel atomic distance (Å)
300	57	59	61
325	64	67	70
350	68	71	75

Table 4

COM angle of coronavirus with a metallic matrix

Temperature (K)	Fe atomic angle	Al atomic angle	Steel atomic angle
300	92	89	93
325	95	88	91
350	88	91	92

Mutual energy of two groups of simulated structures describes their atomic interaction in MD simulations. In this step of our calculations, we report the mutual energy of coronavirus and metallic matrix. Figure 7 shows that the mutual energy between the coronavirus and Fe matrix varies from 32,261/30,569/29,777 to 0 eV by temperature increasing. Zero rates of mutual energy rate show that the coronavirus structure and Fe matrix distance are bigger than the cutoff radius of them at 350 K and these structures departed from each other after 1,000,000-time steps. Coronavirus with Al and steel matrix interacts with a similar manner, and mutual energy of these structures decreases by MD simulation time steps passing (Table 5). Between metallic matrix, steel one has minimum mutual interaction with the coronavirus. That implies the antivirus property of this structure which nominated this metal to antibacterial applications.

Fig. 7

Common energy variation of coronavirus and Fe matrix with MD simulation time steps and temperature

Table 5

Common energy variation of coronavirus and metallic matrix

Temperature (K)	Fe matrix (eV)	Al matrix (eV)	Steel matrix (eV)
300	32,261	28,759	25,222
325	30,569	27,569	22,359
350	29,777	26,335	19,965

Larger simulation time steps, volume of virus structure rises in the presence of steel atoms and so, the atomic distance increases, too (see Fig. 8). More the distance of atoms in the coronavirus structure, the stability of the virus decreases. From Fig. 9 by the time evolution from 0 to 1,000,000-time steps, the volume of coronaviruses increases from 179,091 to 205,283 Å3. Coronavirus volume does not change drastically by simulation time evolution in the presence of Fe and Al matrix (Table 6). Volume of the virus varies from 179,091 to 179,095/179,098 Å3 for Fe/Al matrix, respectively. This atomic behavior indicates that the use of steel in medical applications can have an important role in the elimination of coronavirus.

Fig. 8

Variation of coronavirus volume with time steps: a n = 0, b n = 250,000, n = 500,000, n = 750,000, and n = 1,000,000

Fig. 9

Coronavirus volume variation with metallic matrix and MD simulation time steps

Table 6

Coronavirus volume variation with metallic matrix at T = 300 K

Matrix type	Volume(Å3)
Fe	179,095
Al	179,098
Steel	205,283

Furthermore, with temperature increasing from 300 to 350 K, the final volume of coronavirus varies from 205,283 to 253,244 Å3 in the presence of steel matrix (Table 7). Figure 10 shows the volume variation of the coronavirus structure in the presence of various metallic matrix. From these results, we can say that coronavirus atomic stability do not reduce in the presence of Fe and Al matrix.

Table 7

Coronavirus volume variation with metallic matrix and molecular dynamics simulation temperature

Temperature (K)	Fe matrix (Å3)	Al matrix (Å3)	Steel matrix (Å3)
300	179,095	179,098	205,283
325	179,123	179,269	223,698
350	179,359	181,433	253,244

Fig. 10

Time evolution of coronavirus volume in vicinity of metallic matrix with MD simulation temperature

Conclusion

The present study investigates the atomic interaction between the coronavirus and metallic matrix such as Fe, Al, and steel with MD simulations. In our calculations, the coronavirus is shown by S,O, N, and C atoms. Result of simulations shows that steel matrix has good properties to prevent coronavirus transmission, which can be used for medical purposes. Generally, the results are as following:

Dreiding force field is the appropriate interatomic potential for MD simulation of coronavirus. Numerically, the total energy of our simulated structures converged to − 28,545 eV, − 22,835 eV, and − 17,124 eV after 1 ns at 300 K in the presence of Fe, Al, and steel matrix.

Center of the mass distance of coronavirus and metallic structures increases by the passage of MD simulations time from 45 to 75 Å where the largest rate of this physical parameter occurs for steel matrix.

Generally, by increasing the simulation temperature, the atomic repulsion between the metallic matrix and coronavirus rises.

By increasing the temperature of simulated structures, the angle of the coronavirus does not follow a logical relation

The mutual energy of coronavirus and metallic structures varies from 19,965 to 32,261 eV after 1 ns. The highest rate of this physical parameter calculated for the Fe matrix, which describes maximum interaction between the coronavirus and these metallic structures in our MD simulations.

The volume of the coronavirus increases by MD simulation time passing in the presence of steel matrix. Numerically, coronavirus volume increases from 179,091 to 253,244 after atomic interaction with steel matrix.