Slimane Laref
email: slimane.laref@kaust.edu.sa

Fouzi Harrou
email: fouzi.harrou@kaust.edu.sa

Bin Wang

Ying Sun

Amel Laref

Taous-Meriem Laleg-Kirati

Takashi Gojobori

Xin Gao

Fouzi Harro

Synergy of Small Antiviral Molecules on a Black-Phosphorus Nanocarrier: Machine Learning and Quantum Chemical Simulation Insights

Keywords: MD, machine learning, DFT, inhibitor, black phosphorus, thermodynamic, molecular states, drug vehicles

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The fast development of clinical biomedicine and nanobiotechnology has inspired the current generation of diverse inorganic nanoparticles that provide various pattern modalities as possible substitutions in addressing diseases by synergistic therapy, and principally to treat metastatic pathology [START_REF] Kotagiri | Breaking the depth dependency of phototherapy with Cerenkov radiation and low-radiance-responsive nanophotosensitizers[END_REF][START_REF] Hong | Multifunctional in vivo vascular imaging using near-infrared II fluorescence[END_REF][START_REF] Yang | Nano-graphene in biomedicine: Theranostic applications[END_REF]. Two-dimensional (2D) materials have fascinated extensive consideration due to their exceptional and distinctive properties [START_REF] Ferrari | Cancer nanotechnology: Opportunities and challenges[END_REF][START_REF] Lin | Insights into 2D MXenes for versatile biomedical applications: Current advances and challenges ahead[END_REF][START_REF] Laref | Characterization of graphenefullerene interactions: Insights from density functional theory[END_REF]. Encouraged by the great success of graphene, a variety of 2D materials have been deployed for various property-related applications such as energy, environmental science, catalysis, physics, and biomedicine [START_REF] Backes | Production and processing of graphene and related materials[END_REF][START_REF] Gao | Artificially stacked atomic layers: Toward new van der Waals solids[END_REF][START_REF] Laref | Surface passivation effects on the electronic and optical properties of silicon quantum dots[END_REF][START_REF] Novoselov | A roadmap for graphene[END_REF][START_REF] Wick | Classification framework for graphene-based materials[END_REF][START_REF] Laref | Impact of phosphorous and sulphur substitution on Dirac cone modification and optical behaviors of monolayer graphene for nano-electronic devices[END_REF][START_REF] Smith | properties, and applications of graphene oxide/reduced graphene oxide and their nanocomposites[END_REF][START_REF] Laref | Electronic structure and optical characteristics of AA stacked bilayer graphene: A first principles calculations[END_REF].

For example, transition metal dichalcogenides (TMDs), nitrides and carbonitrides (MXenes), and hexagonal boron nitride (h-BN), and their derivatives are realized with a monolayer and a few layers, which exhibit weak van der Waals interlayer bonding and strong in-plane covalent bonding [START_REF] Huang | Freestanding palladium nanosheets with plasmonic and catalytic properties[END_REF][START_REF] Anasori | 2D metal carbides and nitrides (MXenes) for energy storage[END_REF][START_REF] Laref | A Peculiar Binding Characterization of DNA (RNA) Nucleobases at MoOS-Based Janus Biosensor: Dissimilar Facets Role on Selectivity and Sensitivity[END_REF]. Thereby, the biological behavior of Xenes 2D materials is under outstretched investigation to check their biosafety and biocompatibility, such as biodegradability, cytotoxicity, and toxic derivatives [START_REF] Ma | Degradation-by-design: How chemical functionalization enhances the biodegradability and safety of 2D materials[END_REF]. The broad spectra of potential uses and clinical implementation of the Xenes 2D systems in biomedical drug therapy still require proper control of biodegradability [START_REF] Liu | Semiconducting black phosphorus: Synthesis, transport properties and electronic applications[END_REF] and toxicity, as well as an explicit understanding of the biochemical procedures occurring in the drug molecules on 2D nanomaterials and their interfaces. The first aspect is particularly valuable since the understanding of the orientations of molecules with the multi-layered materials can characterize the process of raising biologically friendly materials for gene silencing, scattering complexes, hybrid pharmaceutic ingredients, and biomedicine sensing. In the present study, motivated by the latest developments in 2D nanomaterial synthesis and their encouraging solicitation in anti-tumor drugs [START_REF] Shen | Fe 3 O 4 nanoparticles in targeted drug/gene delivery systems[END_REF][START_REF] Rahimi | BC 3 graphene-like monolayer as a drug delivery system for nitrosourea anticancer drug: A first-principles perception[END_REF][START_REF] Hashemzadeh | Covalent organic framework as smart and high efficient carrier for anticancer drug delivery: A DFT calculations and molecular dynamics simulation study[END_REF], we conducted a comprehensive investigation in order to use 2D materials as nanocarriers for drug delivery. One such potential pharmacological option for treating various desists is ebselen (EB), while another is favipiravir (FP) [START_REF] Liu | Research and development on therapeutic agents and vaccines for COVID-19 and related human coronavirus diseases[END_REF][START_REF] Jin | Structure of Mpro from SARS-CoV-2 and discovery of its inhibitors[END_REF]. Favipiravir has recently been used to treat infections caused by the Ebola virus, and it was found to be safe during clinical studies [START_REF] Scavone | Current pharmacological treatments for COVID-19: What's next?[END_REF][START_REF] Sreekanth Reddy | Tackling COVID-19 using remdesivir and favipiravir as therapeutic options[END_REF]. During the outbreaks of the SARS-CoV-1 and MERS-CoV infections, favipiravir was once more reviewed, and it was shown to be effective against these viruses in animal model systems [START_REF] Guy | Rapid repurposing of drugs for COVID-19[END_REF][START_REF] Dong | Discovering drugs to treat coronavirus disease 2019 (COVID-19)[END_REF][START_REF] Delang | Favipiravir as a potential countermeasure against neglected and emerging RNA viruses[END_REF][START_REF] Furuta | Favipiravir (T-705), a broad spectrum inhibitor of viral RNA polymerase[END_REF][START_REF] Furuta | T-705), a novel viral RNA polymerase inhibitor[END_REF]. Ebselen is an organoselenium compound with cytoprotective, anti-inflammatory, and anti-oxidant properties. This compound has previously been investigated to treat multiple diseases, including hearing loss [START_REF] Lynch | Development of ebselen, a glutathione peroxidase mimic, for the prevention and treatment of noise-induced hearing loss[END_REF] and bipolar syndrome [START_REF] Singh | A safe lithium mimetic for bipolar disorder[END_REF]. Ebselen has minimal cytotoxicity [START_REF] Renson | 2-Phenyl-1,2-benzisoselenazol-3(2H)-one Containing Pharmaceutical Preparations and Process for the Treatment of Rheumatic Diseases[END_REF], and its safety in humans has been estimated in a number of clinical trials [START_REF] Kil | Safety and efficacy of ebselen for the prevention of noise-induced hearing loss: A randomised, double-blind, placebo-controlled, phase 2 trial[END_REF][START_REF] Masaki | Effects of the potential lithium-mimetic, ebselen, on impulsivity and emotional processing[END_REF][START_REF] Chen | Cinanserin is an inhibitor of the 3C-like proteinase of severe acute respiratory syndrome coronavirus and strongly reduces virus replication in vitro[END_REF][START_REF] Carmo | Clearance and persistence of SARS-CoV-2 RNA in patients with COVID-19[END_REF]. Ebselen has the predominates inhibitor targeting the main protease (Mpro) of SARS-CoV-2 (e.g., Mopar is a crucial enzyme of coronaviruses) [START_REF] Jin | Structure of Mpro from SARS-CoV-2 and discovery of its inhibitors[END_REF]. Recently, in vitro human cell studies reported that both favipiravir and ebselen could effectively stop the replication of SARS-CoV-2 at an earlier stage [START_REF] Jin | Structure of Mpro from SARS-CoV-2 and discovery of its inhibitors[END_REF].

In recent years, the application of machine learning techniques in nanomaterials and pharmaceutical technologies has gained significant attention from researchers. Several studies have used backpropagated neural networks to model various systems [START_REF] Khan | Neural artificial networking for nonlinear Darcy-Forchheimer nanofluidic slip flow[END_REF][START_REF] Shoaib | Intelligent backpropagated neural networks application on Darcy-Forchheimer ferrofluid slip flow system[END_REF], including Darcy-Forchheimer slip flow models for nanomaterials and ferrofluids [START_REF] Raja | Intelligent neuro-computing for entropy generated Darcy-Forchheimer mixed convective fluid flow[END_REF], non-Fourier heat flux models for transient heat exchange [START_REF] Varun Kumar | Backpropagated Neural Network Modeling for the Non-Fourier Thermal Analysis of a Moving Plate[END_REF], and ternary nanofluid flow between parallel plates [START_REF] Sharma | Numerical and Levenberg-Marquardt backpropagation neural networks computation of ternary nanofluid flow across parallel plates with Nield boundary conditions[END_REF]. The accuracy of these models was evaluated using reference datasets obtained via the Homotopy Analysis Method and numerical methods. The effects of non-dimensional parameters on the systems were also analyzed. In this paper, by utilizing accelerated molecular dynamics, machine learning, and quantum chemicals based on first-principles density functional theory (DFT), we briefly report the structural stability and thermochemical properties of FP and EB functionalized within a typical phosphorene 2D monolayer. The chemical bonding of FP and EB in 2D black-phosphorus materials is revealed by considering the vacuum and continuum solvent model. However, by calculating the Gibbs free energy changes of each step, we demonstrated that the temperature and pH values would be used to tune the releasing ratio of the antiviral molecules onto the phosphorene sheet. Our results provide new fundamental insights into exploring innovative 2D materials as new nanocarriers for antiviral drug therapy.

These atomistic findings open the door to the parametrization of multiscale methodologies, such as contemporary machine learning based on electronic theories, continuum models, Monte Carlo simulations, and coarse-grained methods, to address the interplay between the interaction of drugs/2D monolayers with the cell membrane as a function of their morphology. In addition, we rigorously explored the efficiency of four machine learning models to accelerate the prediction of Hamiltonian and non-bond energy from EB, FP, and FP-EB, namely Bagged Trees, Gaussian Process Regression (GPR), Support Vector Regression (SVR), and Regression Trees (RT). The reason for considering multiple models is to compare and evaluate their performance in accurately predicting the behavior of the antiviral molecules in the 2D monolayer. We calibrated the machine learning models using Bayesian optimization based on the training data. Results revealed the superior performance of the GPR model by reaching an averaged R 2 of 0.9649. Training efficient and accurate models for approximating DFT is the final step in using machine learning to gain more insight into the design of new drugs.

The remainder of this paper is organized as follows. Section 2 presents the used data and the obtained results. In Section 3, we discuss the pharmachemical implications. Section 4 briefly describes the preliminary materials, including the Mean Squared Error (MSE), the investigated machine learning models, and ab initio calculations. Finally, we offer conclusions in Section 5.

Outcomes and Discussion

Data Description

Molecular dynamics is a powerful approach that enables atoms in the systems to interact with each other based on given parameters, such as temperature, time, pressure, and other conditions. The evolution of molecular simulations can be presented by various terms in the Hamiltonian (kinetic energy, potential energy, interaction energy, free energy). Thereby, the consistency of MD methods combined with ML and other quantum techniques can be considered. The main key in this work is to improve the Hamiltonian and nonbond energy terms from MD simulations as a starting point. The datasets are generated based on the COMPASS III [START_REF] Akkermans | COMPASS III: Automated fitting workflows and extension to ionic liquids[END_REF] force field of molecular dynamics simulations for small antiviral molecules on BP single layer. By using the ML algorithms mentioned above, the Hamiltonian and the interaction energy are also trained at periodic boundary conditions (PBC) of polycyclic molecules on top of a phosphorene monolayer. Training efficient and accurate models for approximating DFT is the final step in using ML to aid in the design of new drugs.

Data Analysis

This part is dedicated to assessing the capability of the investigated machine learning models in predicting Hamiltonian and non-bond energy from three different datasets, MLFF-EB, MLFF-FP, and MLFF-FP_EB.

Figure 1 displays the distribution of the Hamiltonian energy from three different datasets. Plotting the probability density can be a useful tool for exploratory data analysis. Generally speaking, it can help visualize a data distribution and gain insights into its properties. It shows the relative frequency of different values and the shape of the distribution, which can identify patterns or anomalies. It can also be used to calculate statistical measures such as mean, variance, and skewness, providing further insight into the properties of the data. From Figure 1, we observe that these datasets are non-Gaussian-distributed. It would challenge traditional prediction methods, such as principal component regression, designed based on the Gaussian assumption of data. Thus, machine learning models designed without assumptions about the data distribution could be promising.

Figure 2a-c depicts the pairwise correlation coefficients between the multivariate data of antiviral drugs on the BP monolayer, MLFF-EB, MLFF-FP, and MLFF-FP_EB, respectively. As the calculation proceeds, there is no change in the temperature profile. This suggests that the model is dynamically stable. We observe that the temperature is absolutely correlated with kinetic energy. Once more, the temperature is uncorrelated with the other variables since it is related to the variations of the potential energy, non-bond energy, and total energy, respectively. Furthermore, it is also observed that in the Hamiltonian, potential energy, kinetic energy, non-bond energy, and total energy are constant throughout the simulation processes. We found a lowering in the kinetic energy that leads to general stability. However, this indicates that the model used in this study shows additional thermal stability of the system. The relaxed geometries are further used for creating a black phosphorus sheet for loading antiviral molecules that characterize the drug-BP based nanocarrier. Accordingly, the optimized ML structures (MLFF-EB, MLFF-FP, and MLFF-FP_EB, at T = 0 K) are in good agreement with the interaction energies and van der Waals energies found by means of DFT calculations [START_REF] Laref | Computational Studies of Auto-Active van der Waals interaction Molecules on Ultra-thin Black-Phosphorus Film[END_REF]. For more information, please see Tables S1 andS2, respectively. Finally, during the calculations, the entire cell with its external pressure was fixed, and it had a moderate correlation with respect to bonding energy, angle energy, and TPE, as depicted in Figure 2.

Prediction Results

Hamiltonian Energy Prediction

Three datasets were used to assess the performance of the machine learning models. Each sub-set of data consists of eight input variables (i.e., total kinetic energy, angle energy, torsion energy, inversion energy, van der Waals energy, non-bond energy, pressure, and temperature) and one out-variable (i.e., Hamiltonian energy). Each data set comprises around 2000 samples. In the first experiment, the dataset was divided into two subsets: training composed of 80% and testing with 20%. The k-fold cross-validation technique was considered in constructing these models based on the training data. Specifically, a 5-fold cross-validation procedure was employed in training the investigated models. The Bayesian optimization approach was employed to optimize the GPR, SVR, RT, and BT models to determine the optimal parameters that minimize the Mean Absolute Error between the actual and predicted Hamiltonian energy using the training data.

Table 1 summarizes the evaluation metrics, RMSE, MAE, R 2 , and MAPE, computed from the testing datasets for MLFF-EB, MLFF-FP, and MLFF-FP_EB. Results revealed that the investigated machine learning models exhibited good capability in predicting Hamiltonian energy. We also notice from Table 1 that no single approach uniformly dominates the other models for the three considered datasets. For the MLFF-EB dataset, the SVR model attained the best prediction accuracy with an R 2 of 0.9767. In this case study, the GPR model obtained an accurate prediction with an R 2 of 0.9614 for the MLFF-EB dataset, and the best prediction quality was obtained by the GPR model with an R 2 of 0.9669 when applied to the MLFF-FP_EB dataset.

Table 2 lists the average evaluation metrics per model for Hamiltonian energy prediction. The closest R 2 value to 1 and the lowest RMSE, MAE, and MAPE values characterize the best prediction performance. Results in Table 2 show that all four models performed well in predicting the Hamiltonian energy, with R 2 values ranging from 0.9478 to 0.9649, indicating a good fit between the predicted and actual values. The GPR model had the highest R 2 value (0.9649) and the lowest values for RMSE, MAE, and MAPE, indicating that it performed slightly better than the other models in terms of accuracy. An R 2 value of 0.9649 means that the GPR model can capture 96.49% of the variation in the Hamiltonian energy of the antiviral molecules in the phosphorene monolayer. In addition, the GPR model achieved the lowest RMSE (0.0789) and MAE (0.0521) values among all models, indicating that the model had the smallest average deviation between predicted and actual values. Additionally, the GPR model had the lowest MAPE value (0.0095), indicating that the average percentage error between predicted and actual values was the lowest among all models. Therefore, the GPR model performed the best in terms of all evaluation metrics. This indicates that the model is highly accurate in predicting the Hamiltonian energy and can be used as a reliable tool for designing new antiviral molecules.

Table 1. Prediction quality of the Hamiltonian energy using the four optimized models. It is worth examining the prediction errors from the four investigated models (i.e., BT, GPR, SVR, and RT). The prediction error refers to the deviation between the true Hamiltonian energy measurements and their predicted values. Figure 3a-c displays the boxplots of the prediction errors for each model based on testing datasets from the MLFF-EB, MLFF-FP, and MLFF-FP_EB datasets, respectively. Visually, we observed that the prediction errors of the GPR and SVR models are concentrated around zero, indicating better prediction accuracy than the BT and RT models. Furthermore, it can also be seen that the GPR model with narrower boxes and whiskers reaches superior performance compared to the SVR for the MLFF-FP and MLFF-FP_EB datasets (Figure 3). Moreover, SVR and GPR models show relatively comparable performance under the MLFF-EB dataset (Figure 3).

MLFF-FP

Methods

As a visual illustration, the scatter graphs and plots of actual and predicted Hamiltonian energy of MLFF-EB, MLFF-FP, and MLFF-FP_EB, when applying the trained GPR model to testing data are depicted in Figure 4. To simplify visual readability, we presented only the results from the best model, and the results from the other models were omitted. We observe that the predicted values of Hamiltonian energy obtained by the GPR models are fairly close to the actual data of MLFF-EB, MLFF-FP, and MLFF-FP_EB, which indicates the promising performance of the GPR model.

Non-Bond Energy Prediction

In the second experiment, the main key objective was to predict and improve the non-bond energy for MLFF-EB, MLFF-FP, and MLFF-FP_EB based on seven input variables (i.e., total kinetic energy, angle energy, torsion energy, inversion energy, van der Waals energy, pressure, and temperature) related to ab initio results [START_REF] Laref | Computational Studies of Auto-Active van der Waals interaction Molecules on Ultra-thin Black-Phosphorus Film[END_REF]. Towards this end, we constructed four machine learning models (BT, SVR, GPR, and RT) based on a 5-fold crossvalidation with the training dataset. The parameters of each model have been optimized using the Bayesian optimization approach in training. The trained models have been employed for non-bond energy based on testing datasets. Then, the statistical metrics (R 2 , RMSE, MAE, and MAPE) are computed to compare the prediction performance of the investigated machine learning models (Table 3). For the three datasets, the GPR approach dominates the other approaches (i.e., BT, SVR, and RT) in terms of the four performance evaluation metrics (Table 3). Regarding the MLFF-FP, the four models provide relatively comparable performance with an R 2 around 0.97; the GPR model slightly outperformed the other models with an R 2 of 0.9792. Considering the MLFF-EB data, the GPR model obtained the best prediction performance for MLFF-EB and MLFF-FP_EB with an R 2 of 0.9735 and 0.9725, respectively. Table 4 lists the averaged evaluation metrics per model for non-bond energy prediction. The GPR model performed the best among all the models with an R 2 value of 0.9751, indicating a very good fit between the predicted and actual non-bond energies. The GPR model also had the lowest RMSE and MAE values, indicating that the predicted values were very close to the actual values. The MAPE value for the GPR model was 0.9579%, the lowest among all models, indicating that the model had the lowest percentage error in prediction. Therefore, based on these evaluation metrics, the GPR model was the best performer among all the models for non-bond energy prediction. Figure 5a-c shows the boxplots of the prediction errors for the four models based on testing datasets (MLFF-EB, MLFF-FP, and MLFF-FP_EB). Visually, we observe that the boxplots are around zero, which means that the prediction quality is good and the prediction errors are small. Additionally, it can be seen that the GPR model provides relatively slightly less prediction errors than the other models for the three considered datasets. The actual MD results and the predicted non-bond energy from the four models based on testing datasets are depicted in Figure 6a-c. We found that the predicted values of non-bond energy from the GPR model closely followed the measured data. Figure 6d-f shows the scatter plot of the measured and predicted values of non-bond energy from the GPR model. A scatter plot is a useful tool to visually compare a dataset's predicted and measured data. In Figure 6d-f similarly, it is well-noted that the close link concerning virus and plasma antigen levels of tissue type (plasminogen) has a substantial clinical impact [START_REF] Zuo | Plasma tissue plasminogen activator and plasminogen activator inhibitor-1 in hospitalized COVID-19 patients[END_REF]. Accordingly, as a result of defection tissue, the normal and infected tissue should have different pH numbers. As a result, the stability of different respiratory diseases is the maximum at a pH close to 6, so the application of the voltages also affects the pH of the plasma tissue. The voltage induced is widely used in physical therapy to treat a variety of diseases. However, the pH measurement is simulated in the solvent environment by placing altered charged ions. On the other hand, the solid-solvent interface provides an alternative challenge for biomaterial simulation due to the large number of degrees of freedom hosted by the liquid states. For this purpose, we used an implicit solvation model, where the inhibitor molecule/2D biomaterial is treated at the level of a quantum mechanic state, and the solvent is treated as a continuum model described by the relative permittivity. This implicit solvation model is found to mimic the correct experimental solvation energy in various solvents [START_REF] Cheng | Explanation of dramatic pH-dependence of hydrogen binding on noble metal electrode: Greatly weakened water adsorption at high pH[END_REF].

To predict the pH sensitivity of the binding characteristics, we follow the method in Ref. [START_REF] Cheng | Explanation of dramatic pH-dependence of hydrogen binding on noble metal electrode: Greatly weakened water adsorption at high pH[END_REF] by adding or removing electrons number of valence states, expressed by the modification of the work function and the applied voltage, which has been demonstrated to provide consistent results with the experiments [START_REF] Ou | Plug-and-play nanorization of coarse black phosphorus for targeted chemo-photoimmunotherapy of colorectal cancer[END_REF]. The adsorption energy of neutral FP, EB, and FP + EB on the BP single layer is -2.61, -5.41, and -8.86 meV/Å 2 , respectively; inclusion of 0.1 extra electrons will promote the energy of the system and cause a lowering of the binding energy. For example, adding half electrons in FP, EB, and FP + EB on BP surface changes the energy to 4.02, -1.95, and 2.82 meV/Å 2 , respectively, and the desorption of the hybrid FP + EB drugs becomes spontaneous. This scheme allows us to frequently control the chemical desorption by introducing fractional charge at the interface under the Poisson-Boltzmann approximation. In particular, pH-dependent adsorption depicted in Figure 7 suggests that larger adsorption energy is obtained in a polar environment with a higher pH value. Thus, when they are expected to release pharmaceutical drugs in a pointing position, an application of a low pH value with a smaller positive voltage would be appropriate.

Thermotherapy Properties

The electronic and pH-influencing properties of antiviral drugs were described above.

The key character of chemical adsorption nature related to temperature dependence is the next step in this investigation. Correspondingly, as represented in Figure 8, the temperature indeed alters the adsorption energy of the molecules at the BP monolayer, either at vacuum or solvent level. It can be well-noticed that by considering the vibrational correction, despite the fact that the entropy change is probably dominated by translational entropy rather than vibration, it shifts the adsorption to more positive or less negative numbers at a higher temperature. The zero-point energy at lower temperature values (0-100 K) has a minor effect on the energy (between 1-2 meV/Å 2 for FP, EB, and FP + EB, respectively). On the other hand, increasing temperature from 150 to 420 K causes shifts of up to 4.1 meV/Å 2 for FP, EB, and FP + EB on the BP-vacuum interface with the incorporation of the corrected vibrational free energy. Equivalent trends have been found by adding the continuum solvent level to these antiviral-molecules/BP 2D-system. Surprisingly, the adsorption remains stable up to 450 K, at the vacuum phase for hybrid FP + EB on the BP sheet. The thermotherapy is vigorous, and the releasing rate of the antiviral molecules will be highly enhanced due to the large drop in the adsorption energy at higher temperatures.

In general, the drug will have a larger release rate beyond 350 K, and it might be due to the height's thermal deviation at a larger temperature. In the ex vivo evaluation process, the phosphorene-loaded drugs can be used as an anti-defected cell drug that has a pertinent act in the thermal therapeutic [START_REF] Laref | Computational Studies of Auto-Active van der Waals interaction Molecules on Ultra-thin Black-Phosphorus Film[END_REF]. The temperature could upturn to around 450 K within a few minutes. The excellent therapeutic efficiency of these drugs in the black phosphorus nanocarrier is firmly associated to the large ability to target drug-release and the hyperthermia of the infected cell. Further, these findings suggest that the hybrid FP + EB therapy of low nanocarrier-based delivery systems can be activated effectively [START_REF] Ou | Plug-and-play nanorization of coarse black phosphorus for targeted chemo-photoimmunotherapy of colorectal cancer[END_REF]. This comprehensive investigation will contribute to optimizing virus growth and storage conditions, facilitating the molecular characterization of this important pathogen. The interest of accommodating those inhibitors in close contact with the phosphorene sheet as a new delivery scheme is that it could continuously control the desorption free energy of the aromatic molecules by applying an electric field, together with the presence of external stimuli to tune their pH environment, and changing their temperature.

Materials and Methods

This section is dedicated to molecular dynamics simulations and machine learning models deployed for this study, which is related to describing nanocarrier loaded drugs prediction. ML approaches are effective methods for accelerating DFT relaxations and require significantly fewer iterations (DFT calculations) for convergence. It has shown excellent transferability along FP, EB, and hybrid FP and EB molecules physisorbed on a BP single layer, and closely describes the ab initio interaction energy of the entire system. This is an encouraging initiative toward a common class of force-field and ab initio approaches combined with ML for drug molecules on ultra-thin film crystalline 2D material.

Molecular Dynamics Simulations

Molecular dynamics (MD) simulations based on atomic force fields can accurately simulate the non-covalent geometry of a large number of molecular systems from nonequilibrium systems to thermodynamic states at the atomic level. In the present study, the current data were obtained by means of molecular dynamics simulations, as implemented in Forcite software [START_REF] Akkermans | COMPASS III: Automated fitting workflows and extension to ionic liquids[END_REF]. To properly adsorb FP and EB drugs on a BP single layer, we used a (5 × 4) supercell model of black phosphorus containing 80 phosphorus atoms for the nanocarrier. A vacuum region of 20 Å is introduced perpendicular to the BP film to prevent unphysical periodic effects. A Hamiltonian including kinetic energy, potential energy, interaction energy, and free energy as a function of the change of temperature and pressure was simulated at 350 K and atmospheric pressure, respectively. We used an NVT canonical ensemble based on the Nosé thermostat. During the MD simulations, we adopted a time step of 0.3 fs for an ensemble of 100 ps of calculations to reach the thermodynamic states and accurate MD calculations.

Gaussian Process Regressor

GPR models, which are within nonparametric kernel-driven learning models, showed extended modeling ability for handling nonlinear prediction problems because of their nonlinear approximation capabilities [START_REF] Xie | Gaussian processes for short-term traffic volume forecasting[END_REF][START_REF] Harrou | Monitoring of photovoltaic systems using improved kernel-based learning schemes[END_REF][START_REF] Lee | Wind power prediction using ensemble learning-based models[END_REF]. They are well-known for their ability and flexibility in modeling small-sized data with a Gaussian or non-Gaussian distribution, as well as providing uncertainty measures on predictions [START_REF] Williams | Gaussian processes for regression[END_REF][START_REF] Mackay | Gaussian Processes-A Replacement for Supervised Neural Networks?[END_REF]. Moreover, the GPR model's benefit is its capacity to provide a confidence interval to evaluate the reliability of the predicted values, making it widely used in numerous applications [START_REF] García-Nieto | Prediction of outlet dissolved oxygen in micro-irrigation sand media filters using a Gaussian process regression[END_REF]. The GPR models have been utilized in various applications, including renewable energy systems monitoring [START_REF] Harrou | Monitoring of photovoltaic systems using improved kernel-based learning schemes[END_REF][START_REF] Alkesaiberi | Efficient wind power prediction using machine learning methods: A comparative study[END_REF], COVID-19 spread prediction [START_REF] Alali | A proficient approach to forecast COVID-19 spread via optimized dynamic machine learning models[END_REF], and spatio-temporal PM 2.5 prediction [START_REF] Yang | Prediction of hourly PM 2.5 using a space-time support vector regression model[END_REF].

To introduce GPR, we have to consider the response y of a function f at the input x, which is expressed as [START_REF] Schulz | A tutorial on Gaussian process regression: Modelling, exploring, and exploiting functions[END_REF],

y i = f (x i ) + ε i . (1)
Here, ε ∼ N (0, σ 2 ε ) is an additive noise, and f (x) is considered as a random variable. It is worth pointing out that the uncertainty on f could decrease significantly by the observation of the function's output at different input points.

The function f (x) is considered following a Gaussian process. Accordingly, y i follows a joint Gaussian distribution [START_REF] Schulz | A tutorial on Gaussian process regression: Modelling, exploring, and exploiting functions[END_REF]:

y = [y 1 , y 2 , . . . y n ] ∼ N (m(x), K + σ 2 I), (2)
where m(x) = [m(x 1 ), m(x 2 ), . . . m(x n )] refers to the vector of mean values m(•), I denotes the identity matrix, and K refers to the n × n covariance matrix with (i, j)th element

K ij = k(x i , x j ). Within the GPR framework, K ij = k(x i , x j
) is termed the kernel function [START_REF] Williams | Gaussian processes for regression[END_REF][START_REF] Schulz | A tutorial on Gaussian process regression: Modelling, exploring, and exploiting functions[END_REF]. Numerous kernels have been used in the literature, including the Linear kernel, Rational Quadratic (RQ) kernel, Squared Exponential (SE) kernel, Matern 5/2 (M52) kernel, and Exponential (Exp) kernel. The two commonly used kernels (i.e., linear and squared exponential) are expressed respectively as [START_REF] Williams | Gaussian processes for regression[END_REF]:

k Lin (x i , x j ) = θ 2 1 + θ 2 2 (x i -θ 3 )(x j -θ 3 ). ( 3
)
k SE (x i , x j ) = θ 1 exp x i -x j 2 θ 2 . ( 4
)
To optimize the GPR model in the training stage, we need to determine the kernel parameters that enable maximizing the following likelihood [START_REF] Schulz | A tutorial on Gaussian process regression: Modelling, exploring, and exploiting functions[END_REF][START_REF] Seeger | Gaussian processes for machine learning[END_REF].

θ opt = arg max θ L(θ), (5)
where θ = [θ 1 , θ 2 , . . .] are the parameters of the kernel (also called hyperparameters), the mean values m(.) are taken to be zero, and

L(θ) = 1 (2π) n |K + σ 2 I| exp - 1 2 (y (K + σ 2 I)y) . (6)
Many efforts have been reported in the literature to calibrate the hyper-parameters by using grid search or random search. Here, Bayesian optimization will be adopted to fine-tune the GPR hyperparameters by maximizing the marginal likelihood in [START_REF] Lin | Insights into 2D MXenes for versatile biomedical applications: Current advances and challenges ahead[END_REF] with respect to θ [START_REF] Nguyen | Applying Bayesian Optimization for Machine Learning Models in Predicting the Surface Roughness in Single-Point Diamond Turning Polycarbonate[END_REF].

We assume x * is a new test input, then the estimated mean and variance associated with ŷ * = f (x * ) = f * are [START_REF] Schulz | A tutorial on Gaussian process regression: Modelling, exploring, and exploiting functions[END_REF][START_REF] Seeger | Gaussian processes for machine learning[END_REF]:

ŷ * = k * (K + σ 2 I) -1 y, (7)
and

Σ * = k * * -k * (K + σ 2 I) -1 k * . (8)
respectively. Then, y * follows a conditional distribution which has the form:

y * |y ∼ N ( ŷ * , Σ * ), (9)
where K = k(X, X), K * * = k(X * , X * ), and K * = k(X, X * ) are the covariance matrices computed based on the training data, the testing set, and both training and test sets, respectively. Overall, GPR-based prediction involves two main steps. The first step is to estimate the hyperparameters, which include the kernel and noise variance, in the GPR model by maximizing the negative log marginalized likelihood based on the training data. The second step is to compute the predictive posterior distribution of GPR for new inputs, X * , by using the mean of the predictive distribution in Equation ( 8) as the predicted value. For more details about the GPR model, see Refs. [START_REF] Xie | Gaussian processes for short-term traffic volume forecasting[END_REF][START_REF] Harrou | Monitoring of photovoltaic systems using improved kernel-based learning schemes[END_REF][START_REF] Schulz | A tutorial on Gaussian process regression: Modelling, exploring, and exploiting functions[END_REF][START_REF] Murphy | Machine Learning: A Probabilistic Perspective[END_REF].

Bayesian Optimization Procedure

As discussed before, we need to calibrate hyperparameters during training to obtain good prediction performance when constructing a GPR model. Fine-tuning a machine learning model entails modifying its hyperparameters to optimize performance on a particular task. Optimal hyperparameter configuration has a direct impact on the model's performance. Grid Search, Random Search, and Bayesian Optimization are popular methods for hyperparameter tuning in machine learning. Grid Search is a simple but computationally expensive method that exhaustively searches all possible combinations of hyperparameters within a predefined range. While Grid Search can guarantee to find the optimal combination of hyperparameters, it can be very time-consuming and computationally expensive, especially when the number of hyperparameters and their possible values are high [START_REF] Bergstra | Random search for hyper-parameter optimization[END_REF]. Random Search is another method for hyperparameter tuning that randomly samples a set of hyperparameters from a predefined range. Unlike Grid Search, Random Search does not guarantee finding the optimal combination of hyperparameters, but it is generally faster and more computationally efficient than Grid Search [START_REF] Bergstra | Random search for hyper-parameter optimization[END_REF]. On the other hand, Bayesian Optimization (BO) is a more advanced method that uses probabilistic models and previous evaluations to guide the Search [START_REF] Protopapadakis | An investigation on multi-objective optimization of feedforward neural network topology[END_REF]. It can quickly find good combinations of hyperparameters with fewer evaluations than Grid Search and Random Search.

This study adopts Bayesian optimization to find the optimal values of GPR's hyperparameters [START_REF] Protopapadakis | An investigation on multi-objective optimization of feedforward neural network topology[END_REF]. Recently, the BO approach has been frequently employed to fine-tune hyperparameters in machine learning methods. This could be attributed to its design using Gaussian processes and Bayesian inference to reach global optimization [START_REF] Nguyen | Applying Bayesian Optimization for Machine Learning Models in Predicting the Surface Roughness in Single-Point Diamond Turning Polycarbonate[END_REF]. Importantly, the BO procedure considers the previous evaluations to determine the hyperparameter set and then evaluates the next step, which enables reducing the time of the optimization procedure [START_REF] Shahriari | Taking the human out of the loop: A review of Bayesian optimization[END_REF]. Moreover, the BO procedure can optimize functions without closed form [START_REF] Snoek | Practical bayesian optimization of machine learning algorithms[END_REF]. The optimization process based on the BO procedure needs fewer iterations compared to a grid search optimization.

The main idea of the BO procedure consists of building a probabilistic proxy model for the objective function, employing prior experiments' results as training data [START_REF] Cui | A novel framework of graph Bayesian optimization and its applications to real-world network analysis[END_REF][START_REF] Nikolaidis | Gaussian process-based Bayesian optimization for data-driven unit commitment[END_REF]. Essentially, the proxy model (e.g., the Gaussian process) is not expensive to calculate and enables getting pertinent information on where we should evaluate the true loss function to reach suitable results. Eventually, we consider now the problem of adjusting m hyperparameters P = p 1 , . . . , p m . To this end, the aim is to determine

P * = arg min P g(P|{(x i , y i )} n i=1 ), ( 10
)
where g denotes the the objective function [START_REF] Alali | A proficient approach to forecast COVID-19 spread via optimized dynamic machine learning models[END_REF]. The optimization procedure is controlled by an appropriate acquisition function that decides the next set of hyperparameters to be evaluated [START_REF] Springenberg | Bayesian optimization with robust Bayesian neural networks[END_REF]. Figure 9 depicts the basic idea of the BO procedure to optimize the GPR model during the training phase. Specifically, the mean squared error (MSE) between the actual molecular dynamics data and the GPR predictions is computed at each iteration. The optimized model is obtained once the MSE converges to a small value, close to zero.

SVR Models

Here, we briefly introduce another widely used approach, the SVR model, a flexible data-based approach with good learning capacity via kernel tricks. The key idea underlying the SVR consists in mapping the train data to a higher dimensional space and conducting linear regression in that space. This SVR model is known to be efficient in dealing with nonlinear regression by using the kernel trick, which enables mapping the input features into high-dimensional feature spaces [START_REF] Yu | Support vector regression for real-time flood stage forecasting[END_REF][START_REF] Hong | SVR with hybrid chaotic genetic algorithms for tourism demand forecasting[END_REF]. Moreover, the relevant concept used in designing the SVR model lies in structural risk minimization. It is demonstrated that SVR provides satisfactory performance with limited samples [START_REF] Smola | A tutorial on support vector regression[END_REF]. Thus, VVR models have been broadly exploited in numerous applications, such as solar irradiance prediction [START_REF] Lee | Reliable solar irradiance prediction using ensemble learning-based models: A comparative study[END_REF], wind power prediction [START_REF] Lee | Wind power prediction using ensemble learning-based models[END_REF], and anomaly detection [START_REF] Harrou | Monitoring of photovoltaic systems using improved kernel-based learning schemes[END_REF]. This study will use optimized SVR via Bayesian optimization for comparison.

Bagged Tree Model

Likewise, we discuss another important machine learning model called the bagged tree (BT), also known as bootstrap aggregating [START_REF] Zhang | A gradient boosting method to improve travel time prediction[END_REF]. The essence of this model is based on merging the benefit of the bagging technique and decision trees to enhance prediction quality. Importantly, in the BT model, bootstrap sampling is employed to generate numerous samples from the original dataset. Next, multiple distinct decision trees are built, and their outputs are aggregated to obtain the final output [START_REF] Bauer | An empirical comparison of voting classification algorithms: Bagging, boosting, and variants[END_REF]. Hence, the prediction quality of the decision trees will be improved, the error will be reduced, and the overfitting issue in individual trees will be bypassed significantly [START_REF] Harrou | Wind power prediction using bootstrap aggregating trees approach to enabling sustainable wind power integration in a smart grid[END_REF][START_REF] Ribeiro | Ensemble approach based on bagging, boosting and stacking for short-term prediction in agribusiness time series[END_REF].

Evaluation Metrics

In this study, we assess the accuracy of the predicting models using three metrics: root mean square error (RMSE), mean absolute error (MAE), mean absolute percentage error (MAPE), and R 2 .

• RMSE: Root Mean Squared Error measures the differences between predicted and actual values. It is calculated by taking the square root of the average of the squared differences between the predicted and actual values.

RMSE = 1 n n ∑ t=1 (y t -ŷt ) 2 , (11)
• MAPE: Mean Absolute Percentage Error is a measure of the accuracy of a model in predicting values. It is calculated as the average of the absolute percentage difference between predicted and actual values.

MAPE = 100 n n ∑ t=1 y t -ŷt y t %, (12)
• MAE: Mean Absolute Error is a measure of the average magnitude of the errors in a set of predictions, without considering their direction. It is calculated as the average of the absolute differences between predicted and actual values.

MAE = ∑ n t=1 |y t -ŷt | n , (13)
• R-squared: R 2 is a statistical measure that represents the proportion of the variance for a dependent variable that is explained by an independent variable or variables in a regression model. It is also known as the coefficient of determination, and its value ranges between 0 and 1.

R 2 = ∑ n t=1 [(y t -ȳ) • ( ŷt -ȳ)] 2 ∑ n t=1 (y t -ȳ) 2 • ∑ n t=1 ( ŷt -ȳ) 2 , ( 14
)
where y t is the number of Hamiltonian data, ŷt is its corresponding forecasted Hamiltonian and non-bond energy, and n is the number of records. Importantly, a higher R 2 value indicates a better fit of the model to the data. For example, if the R 2 value is between 0.7 and 0.9, it is generally considered to indicate a good fit. An R 2 higher than 0.9 indicates a very good fit of the model to the data. It means that the model can explain a large proportion of the variability in the data. on the other hand, a lower MAPE value indicates better model prediction accuracy. Lower RMSE, MAE, and MAPE values would imply better precision and prediction quality. Generally, a MAPE value of less than 10% is considered good, while a value greater than 10% is considered poor. However, the threshold for an acceptable MAPE value can vary depending on the application and the specific domain. Therefore, it is important to consider the context and specific requirements when interpreting the MAPE value.

Prediction Framework

This work investigates the feasibility of machine learning methods to predict the Hamiltonian and non-bond energy terms from MD simulations. The prediction methodology employed in this study is depicted in Figure 10. We compared four machine learning methods, including BT, SVR, GPR, and RT. The prediction framework consists of two key steps: model construction and model prediction. The data were first preprocessed to a zero mean and unit variance by subtracting the mean, µ, and dividing by the standard deviation, σ.

y s = (y -µ)/σ. ( 15
)
After obtaining the predicted values, we reversed this step. The normalized data were then split into training and testing sets. The training set was used to construct the machine-learning models. Here, we employ Bayesian optimization to determine the optimal parameters of the machine learning methods. We applied a five-fold cross-validation technique in training the investigated models. After that, the previously constructed models were then utilized to predict the Hamiltonian (or non-bond energy) in the testing step.

Ab Initio Calculations

One of the most popular approaches in quantum chemistry is density functional theory (DFT) that computes a wide variety of properties of almost any kind of atomistic system, such as molecules, crystals, surfaces, interfaces, and even biomaterial carriers when combined with machine learning. Thus, as the next step of this study, all DFT calculations were carried out using the Vienna Ab Initio Simulation Package (VASP 5.4) [START_REF] Kresse | Efficiency of ab initio total energy calculations for metals and semiconductors using a plane-wave basis set[END_REF][START_REF] Kresse | Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set[END_REF]. All electrons with a projected augmented wave (PAW) formalism were used to model the electron-ion interactions [START_REF] Kresse | From ultrasoft pseudopotentials to the projector augmented-wave method[END_REF]. The exchange-correlation functional contribution to the total energy was modeled using the generalized gradient approximation (GGA). In addition to the modified Perdew-Burke-Ernzerhof (PBE) [START_REF] Kresse | From ultrasoft pseudopotentials to the projector augmented-wave method[END_REF], for treating the van der Waals (vdW), interaction corrections through the D3-BJ approach were incorporated in all calculations [START_REF] Perdew | Generalized gradient approximation made simple[END_REF][START_REF] Grimme | A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu[END_REF], as established by Grimme et al. [START_REF] Grimme | Effect of the damping function in dispersion corrected density functional theory[END_REF][START_REF] Grimme | Accurate description of van der Waals complexes by density functional theory including empirical corrections[END_REF]. We used a (5 × 4) supercell model of a black phosphorus monolayer containing 80 phosphorus atoms for the nanocarrier to adsorb FP and EB drugs. A vacuum region of 20 Å was introduced perpendicular to the BP film to prevent unphysical periodic effects. Given the large size of the supercell, numerical integration was employed over the Brillouin zone using a 3 × 3 × 1 Γ-centered k-point grid [START_REF] Grimme | Semiempirical GGA-type density functional constructed with a long-range dispersion correction[END_REF]. A denser 5 × 5 × 1 mesh was applied during the self-consistent field (SCF) for electronic and frequency calculations. The pseudo-wave functions were expanded in a plane-wave basis set with an energy cut-off of 450 eV. The accuracy was ensured by adopting energy convergence criteria of 1 × 10 -5 eV. Therefore, structural relaxation was performed using the conjugate gradient technique until the force on each atom was below 0.01 eV/Å. This approach has been successfully employed in the past to describe fullerene adsorption on single-layer graphene [START_REF] Laref | Characterization of graphenefullerene interactions: Insights from density functional theory[END_REF].

Conclusions

We have systematically unveiled the interaction characteristics of inhibitor drugs on the surface of a phosphorene sheet by using accelerated molecular dynamics simulations, machine learning, and density functional theory calculations. Importantly, four machine learning models (i.e., BT, GPR, SVR, and RT) were employed to predict Hamiltonian and non-bond energy in EB, FP, and FP + EB. Results indicate the promising prediction performance of machine learning models. Specifically, the optimized GPR model obtained the best performance in this case study compared to the other investigated models. Specifically, the GPR model consistently outperformed the other three models (BT, SVR, and RT) for both Hamiltonian and Non-Bond energy predictions. In particular, the GPR model achieved the highest R 2 value of 0.9649 for Hamiltonian energy prediction and 0.9751 for Non-Bond energy prediction. Additionally, the GPR model had the lowest RMSE, MAE, and MAPE values for both energy predictions, indicating its high accuracy and precision in predicting the energy values. Therefore, the GPR model can be considered the best-performing model among the tested machine-learning models in this study. Further, regarding DFT calculations, we carried these out in a vacuum and an aqueous phase to evaluate the binding ability and thermochemical properties for active drug delivery. We found that the hybrid drugs are physisorbed and enable functionalization of the 2D phosphorene, and they show robust thermostability. All the drug molecules displayed strong van der Waals interactions when combined with the phosphorene sheet. Moreover, the calculated pH values and Gibbs free energy changes at different temperatures suggest that the FP and EB enabled the release of the drug from the 2D monolayer at high temperatures and different pH conditions. These findings would point to drug/2D biosystems as a new and potentially effective way to steer drug vehicle delivery in thermotherapy.

Supplementary Materials:

The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/molecules28083521/s1, Figure S1

Figure 1 .

 1 Figure 1. Distribution of Hamiltonian energy: (a) MLFF-EB, (b) MLFF-FP, and (c) MLFF-FP_EB.

Figure 2 .

 2 Figure 2. Cont.

Figure 2 .

 2 Figure 2. Heatmap of the correlation matrix: (a) MLFF-EB, (b) MLFF-FP, and (c) MLFF-FP_EB.

Figure 3 .

 3 Figure 3. Prediction Distribution of prediction errors of the four models based on testing datasets: (a) MLFF-EB, (b) MLFF-FP, and (c) MLFF-FP_EB.

Figure 4 .

 4 Figure 4. Predicted Hamiltonian energy from the GPR model based on test data: (Top) MLFF-EB, (Middle) MLFF-FP, and (Bottom) MLFF-FP_EB.

Figure 5 .

 5 Figure 5. Prediction errors of the prediction models from the testing datasets: (a) MLFF-EB, (b) MLFF-FP, and (c) MLFF-FP_EB.

  , the colored dots represent the data points, where the x-axis represents the predicted values and the y-axis represents the measured non-bond energy values. If the dots are closely clustered around the diagonal line, the model performs well and makes accurate predictions. Thus, from Figure 6d-f, we observe that the predictions are close to the measured values.

Figure 6 .

 6 Figure 6. Predicted non-bond energy from the GPR model based on test data.3. Pharmachemical Implications3.1. Drug Release 3.1.1. The pH Sensitivity The 2D biomaterials are extensively applied in thermotherapy by tuning the drugvirus interactions. It is important to assign an effective nanocarrier for drug therapy;

Figure 7 .

 7 Figure 7. The calculated adsorption energy process as a function of net charge at the solvent phase of the BP sheet with favipiravir, ebselen, and favipiravir + ebselen hybrid drugs. The minus red and plus blue sign colors represent the difference between pH < 0 and pH > 0 values, respectively.

Figure 8 .

 8 Figure 8. The free Gibbs adsorption energy process as a function of temperature at the vacuum and solvent phases of the BP sheet with favipiravir, ebselen, and favipiravir+ebselen hybrid drugs.

Figure 9 .

 9 Figure 9. Schematic drawing of the main idea of the BO-based optimized GPR strategy.

  The model's accuracy was checked by comparing measured data to predicted data via the score indicators: R 2 , RMSE, MAE, and MAPE. An R 2 closer to 1 and lower RMSE, MAE, and MAPE values reflect an accurate prediction.

Figure 10 .

 10 Figure 10. Illustration of the used prediction framework and general steps in antiviral drug discovery for small molecules on black phosphorus, developing quantitative structure-activity relationship models to optimize lead structures for pharmaceutical implication.

  : Molecular dynamics simulations potential energy prediction: (a) snapshot of FP-EB adsorption molecule on Black-Phosphorus, and (b) potential energy change of FP-EB on BP surface. Table S1: ML, and PBE-D3(BJ) non-bond energy prediction of MLFF-EB, MLFF-FP, and MLFF-FP_EB adsorption states in (eV) and (kcal/mol) at

Table 2 .

 2 Averaged evaluation metrics per model for Hamiltonian energy prediction.

		R 2	RMSE	MAE	MAPE
	BT	0.9743	0.0266	0.0217	0.0049
	GPR	0.9663	0.0304	0.0238	0.0053
	SVR	0.9767	0.0253	0.0208	0.0047
	RT	0.9674	0.0299	0.0238	0.0053
			MLFF-EB
	Methods	R 2	RMSE	MAE	MAPE
	BT	0.9513	0.1001	0.0655	0.0115
	GPR	0.9614	0.0891	0.0606	0.0106
	SVR	0.9447	0.1067	0.0809	0.0142
	RT	0.9501	0.1014	0.0654	0.0114
			MLFF-FP_EB
	Methods	R 2	RMSE	MAE	MAPE
	BT	0.9489	0.1459	0.1011	0.0179
	GPR	0.9669	0.1173	0.0718	0.0127
	SVR	0.9665	0.1182	0.0710	0.0125
	RT	0.9259	0.1756	0.1109	0.0196

Table 3 .

 3 Prediction quality of the non-bond energy.

			MLFF-FP
		R 2	RMSE	MAE	MAPE
	BT	0.9772	0.0279	0.0218	0.2466
	GPR	0.9792	0.0266	0.0208	0.2348
	SVR	0.9772	0.0279	0.0211	0.2382
	RT	0.9775	0.0277	0.0221	0.2502
			MLFF-EB
		R 2	RMSE	MAE	MAPE
	BT	0.9685	0.0758	0.0493	0.4496
	GPR	0.9735	0.0695	0.0468	0.4277
	SVR	0.9705	0.0734	0.0492	0.4503
	RT	0.9729	0.0703	0.0486	0.4448
			MLFF-FP_EB
		R 2	RMSE	MAE	MAPE
	BT	0.9648	0.1253	0.0843	2.5297
	GPR	0.9725	0.1107	0.0732	2.2111
	SVR	0.9710	0.1138	0.0771	2.3273
	RT	0.9689	0.1178	0.0795	2.3512

Table 4 .

 4 Averaged evaluation metrics per model for non-bond energy prediction.

		R 2	RMSE	MAE	MAPE
	BT	0.9701	0.0763	0.0518	1.0753
	GPR	0.9751	0.0689	0.0469	0.9579
	SVR	0.9729	0.0717	0.0491	1.0053
	RT	0.9731	0.0719	0.0501	1.0154

T = 0 K. Table S2: ML, and PBE-D3(BJ) van der Waals energy prediction of MLFF-EB, MLFF-FP, and MLFF-FP_EB adsorption states in (eV) and (kcal/mol) at T = 0 K.

Funding:

The authors were supported by King Abdullah University of Science and Technology (KAUST) through Award No. FCC/1/1976-09-01 from the Office of Sponsored Research (OSR). For computer time, this research used the HPC resources of the Supercomputing Laboratory at KAUST.

Data Availability Statement:

The data and results associated with this article are available at https: //github.com/SLIM23-CBRC/MD-ML-Hamiltonian (accessed on 1 April 2023).

Conflicts of Interest:

The authors declare no conflict of interest.


Axelle Moreau

Amelie Bethencourt

Valentin Payet

Mégane Turina

Julia Moulinard

Henri Chabrol

Emeline Chauchard
email: emeline.chauchard@univ-nantes.fr

Amelie Bethencourt

Payet B Valentin

Turina B Julia Mégane

Moulinard

Henri Chabrol

Rage in Video Gaming, Characteristics of Loss of Control Among Gamers

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

ABSTRACT

Online forums devoted to video games dedicate discussions to a phenomenon called "rage" that produces difficulties for gamers. Rage is described as expressing violent and impulsive reactions through excessive frustration or bad faith while gaming. The aim of this qualitative study was to describe rage as it occurs in online video games and to identify its etiology and consequences. A qualitative study included 24 gamers aged 18 to 30 years. Semi-directive interviews were conducted and participants completed an online survey. Two researchers rated each interview using the thematic analysis method and reported an agreement rate of 83%. The gamers described rage as a change in the manner of playing caused by a feeling of frustration during the game. Rage has components that are behavioral (aggressive behaviors), cognitive (poor decision-making) and emotional (anger). Rage can be contagious among gamers and may result in the loss of the game and a "rage quit", an abrupt cessation of the game. Rage appears of prime importance to regulate rage episodes and targeting emotion regulation which can be used as a hook for the prevention of the problematic gaming.

Public Policy Relevance Statement

Rage is a phenomenon of behavioral, emotional and cognitive loss of control that occurs while playing an online video game. Rage could have consequences on psychological wellbeing of players and their teammates. Players seek to control rage episodes to perform in the game; rage could also be a gateway for prevention

Rage in video gaming, characteristics of loss of control among gamers: A qualitative study

In the last two decades, video games have become increasingly included in our everyday life through various technologies (game consoles, laptops, smartphones, or tablets).

In the USA, video games are used in 67% of homes; about 60% of video gamers are under 35 years old and 40% of gamers are female (Entertainment Software Association., 2017). In this age group, gamers mostly play on computers (Entertainment Software Association., 2017; [START_REF] Hakim | Le marché du jeu vidéo en 2014[END_REF]. In France, 73% of the population plays video games and 49% of them play online [START_REF] Hakim | Le marché du jeu vidéo en 2014[END_REF]. Despite the massive use of video games in our societies, few studies have focused on gamers' behaviors and experiences during gaming [START_REF] Kuss | Internet gaming addiction : A systematic review of empirical research[END_REF].

Most of the players have a recreational use of video games, but some individuals may develop excessive use (of over 30 hours a week) or an addictive behavior with consequences on health (American Psychiatric Association, 2013; [START_REF] Batthyány | Computer game playing : Clinical characteristics of dependence and abuse among adolescents[END_REF][START_REF] Porter | Recognizing problem video game use[END_REF][START_REF] Rehbein | Prevalence and risk factors of video game dependency in adolescence : Results of a German nationwide survey[END_REF]. There is currently a debate on gaming disorder, the DSM-5 does not include it as a disorder, whereas the World Health Organization (WHO), in the ICD-11, does recognize it as such . A gaming disorder is characterized by impaired control over gaming, increasing priority given to gaming, and the continuation or escalation of gaming despite the occurrence of negative consequences and functioning impairments over a 12-month period (World Health Assembly, 2018). Gaming disorder risk factors have been identified; some are related to the structure of the game and others to the gamer's vulnerability [START_REF] Caplan | Problematic Internet use and psychosocial wellbeing among MMO players[END_REF][START_REF] Kim | The relationship between online game addiction and aggression, self-control and narcissistic personality traits[END_REF][START_REF] Kuss | Internet Addiction : A Systematic Review of Epidemiological Research for the Last Decade[END_REF][START_REF] Kuss | Internet gaming addiction : A systematic review of empirical research[END_REF][START_REF] Rehbein | Prevalence and risk factors of video game dependency in adolescence : Results of a German nationwide survey[END_REF]. Thus, playing online and playing Massively Multiplayer Online Role-Playing Game (MMORPG), having a strong desire to improve gaming skills, having a powerful identification with the avatar, or a low self-esteem, are considered to increase the risk of developing a gaming disorder or excessive use of video games. Similarly, emotion regulation is involved in the vulnerability of gaming disorders and in negative outcomes of gaming [START_REF] Hemenover | Video games, emotion, and emotion regulation : Expanding the scope[END_REF][START_REF] Marchica | The Role of Emotion Regulation in Video Gaming and Gambling Disorder : A Systematic Review[END_REF].

Gamers need to regulate unpleasant emotions, such as frustration and anger. They are unpleasant feelings that can lead to aggressive behavior. In the framework of video games, anger is mostly studied in the context of violent video games and is associated with aggressive behavior. The assumption that violent video games increase anger and generate aggression yields contrasting results and is and there is currently no consensus about the nature of their association [START_REF] Lee | Violent Video Games and Aggression : Stimulation or Catharsis or Both? Cyberpsychology[END_REF][START_REF] Morrison | Sex differences in anger states after violent video game play and the relationship of empathy and emotional invalidation as preaggressive predictors in young adults[END_REF]. However, anger and anger related ruminations are involved in the relationship between violent video games and aggressive behavior [START_REF] Morrison | Sex differences in anger states after violent video game play and the relationship of empathy and emotional invalidation as preaggressive predictors in young adults[END_REF][START_REF] Yao | Violent video games exposure and aggression : The role of moral disengagement, anger, hostility, and disinhibition[END_REF][START_REF] Zheng | Online violent video games and online aggressive behavior among Chinese college students : The role of anger rumination and self-control[END_REF]. Although anger is mainly discussed in the literature in relation with aggressive behavior and violent video games, anger and rage also appear in nonviolent video games and on forums and platforms of discussion. The mainstream gamers report the experience of violent, impulsive reactions, and loss of control during the game. While it is documented by gamers and described as a part of the game experience there is a lack of literature on these episodes of loss of control during the game, whether it is violent or not (Wu et al., 2021).

An observation of forums dedicated to online video games found that players often mention a particular experience that causes difficulties, named 'rage' or 'tilt'. Rage is defined by players as the expression of violent and impulsive reactions through excessive frustration or bad faith while gaming. Rage appears to cause a change in the gamer's behavior, a loss of control over behavior and emotion, and a loss or a cessation of the game named "rage-quit" [START_REF] Przybylski | Competence-impeding electronic games and players' aggressive feelings, thoughts, and behaviors[END_REF]. Recently, two qualitative studies explored the characteristics of rage among children and adolescents engaged in e-sport [START_REF] Kahila | If the game does not work, it is lagging, or you die in game, you just get furious[END_REF]Wu et al., 2021). The authors described rage as a negative emotional reaction, mainly frustration and anger, associated to a cognitive confusion and a deviance of the usual game play (Wu et al., 2021). To regulate rage, gamers start by quitting the game either to calm down, or to express emotions by crying, shouting, and breaking the materiel, or finally by doing nothing and accepting the rage episode. The main triggers of rage are other players (friend, member of team, or a stranger), oneself by judging his own level of play as low (low self-efficacy and failure) and technical problem that occurred during the game (with material, game, or Internet connection) [START_REF] Kahila | If the game does not work, it is lagging, or you die in game, you just get furious[END_REF]Wu et al., 2021). In these samples of children and adolescents in the context of e-sport, the authors mainly explored the reasons, the emotions and their regulation, as well as the behavior associated with rage. Noteworthy, cognitive functioning and consequences of rage remain unexplored.

Based on the descriptions identified on the forum, rage in video games is experienced by most regular gamers. To our knowledge, no studies focus on describing video game rage among young adults. We based our study on a sample of young adults aged between 18-35 years old given that it is the age group which gather the most players while looking at online data about video games (Entertainment Software Association, 2017), and [START_REF] Bilgihan | A quantitative study exploring the difference between gaming genre preferences[END_REF] argue for an age limit of 35 years old, that separates players according to time spend playing (younger play more) [START_REF] Bilgihan | A quantitative study exploring the difference between gaming genre preferences[END_REF]. This study is exploratory and based on the Interaction of Person-Affect-Cognition-Execution model for addictive behavior [START_REF] Brand | Integrating psychological and neurobiological considerations regarding the development and maintenance of specific Internet-use disorders : An Interaction of Person-Affect-Cognition-Execution (I-PACE) model[END_REF][START_REF] Brand | The Interaction of Person-Affect-Cognition-Execution (I-PACE) model for addictive behaviors : Update, generalization to addictive behaviors beyond internet-use disorders, and specification of the process character of addictive behaviors[END_REF], in order to understand the cognitive, behavioral and emotional characteristics of rage [START_REF] Przybylski | Competence-impeding electronic games and players' aggressive feelings, thoughts, and behaviors[END_REF]. We focus on participants' description of their own rage episodes and their acts, feelings and thoughts, before, during and after rage episodes. The first aim of this exploratory study was to define those episodes of behavioral loss of control, to identify their characteristics and to collect the cognitive, behavioral and emotional data from the gamer's experience. The second objective was to identify the strategies implemented by the players to regulate those episodes.

Methods

Participants

To participate, gamers had to be aged between 18 to 35 years old and had to be online gamers on computers, as well as MMORPG, First-Person Shooter (FPS) or Multiplayer Online Battle Arena (MOBA) players, for at least one hour a week. Twenty-four interviews were conducted, the average length was 12 minutes, fluctuating between 7 to 20 minutes.

Most of the participants were male (n=20).

Measures

Semi-structured interviews

A semi-structured interview guide was created in order to collect data on the experience of raging, including its behavioral, cognitive, emotional and environmental characteristics. It was based on the semi-structured interview guide used in the qualitative research on online poker tilt [START_REF] Moreau | Le TILT au Poker en ligne : Un comportement pathologique transitoire?[END_REF]. The semi-structured interview guide is available in the Supplementary Material.

Thematic analysis

To analyze the interviews, we used a thematic analysis procedure; four Master's students participated in the analyzes. First, each individual student carried out a preliminary thematic analysis of three interviews. Secondly, they pooled the matrix they had made after the preliminary thematic analyzes to combine the themes identified by each assessor and build the rating matrix together. Thirdly, the four investigators rated the same interview individually using the rating matrix to verify and adjust the evaluations. Lastly, the thematic analyzes of the 24 interviews were performed; two students worked on the evaluations of each interview, and each of them rated 12 interviews. The inter-rater reliability was verified and showed a good agreement rate, as 83% of the occurrences were classified identically by the two pairs of students. An occurrence corresponds to a verbatim that fit with one of the thematic identified in the rating matrix.

Online questionnaire

Participants completed an online questionnaire that included sociodemographic and video game use questions, as well as the short version of the Game Addiction Scale (GAS) [START_REF] Lemmens | Development and Validation of a Game Addiction Scale for Adolescents[END_REF]. The GAS evaluates video game addiction with seven items measuring strength, tolerance, humor modifications, relapse, avoidance, conflict and problems that the subject experiences in relation to gaming. Participants indicated the frequency of the behaviors described on a five-point Likert scale (from 1 "never" to 5 "very often"). The scores ranged between 7 to 35, and Cronbach's a was =.82. Use is considered problematic when participants obtain a score of above 3 (sometimes) for at least four items [START_REF] Lemmens | Development and Validation of a Game Addiction Scale for Adolescents[END_REF].

Procedure

Data were collected between December 2016 and March 2017. We called for participants by displaying the study announcement on online social networks and in various places or events (Game Show, stores, University). The individuals interested in participating contacted the research team by email to agree on an appointment. We asked them to complete an online informed consent form. For convenience reasons, the interviews were made by an investigator using Skype. After the interview, the participant was emailed the link to complete the online questionnaire. The investigator was available by chat during the time the participants completed the online questionnaire.

Results

Descriptive statistics

The sample was mostly male and the average age was 22 years old (SD=3.3). 12% (n=3) played once a day; 45% (n=11) played several times a day. The duration of the gaming sessions was more than 2 hours (75%, n=18). The types of games most played by the participants were MOBA (42 %, n=10) and FPS (46%, n=11). The problematic use of video games involved 66.5% (n=16) of the participants. The descriptive data of the sample is presented in Table 1.

Description of the rage phenomenon

All of the participants (N=24) indicated experiencing changes in behavior during video gaming. The terms most commonly used to describe this phenomenon are "rage" (n=24) and "tilt" (n=14). These episodes seem to impact the gamers on an emotional, cognitive and behavioral level. This phenomenon is generally described as temporary loss of control of playing behavior, linked to anger and frustration. "It's a loss of control, in the heat of the moment, that leads to loss of self-control, which is linked to certain events in the game, to a lack of skill, etc."; "Raging is more like getting angry with the other people, being more negative and sometimes even leaving in the middle of a game round"; "I'll get angry at my PC, I'll probably shout out insults (laughs). Erm, shout a lot, even. Then again, I'm not personally the type to throw my controller around, well, maybe a bit in fact."; "The name speaks for itself, raging is really when you blow a fuse. It's when you can't control yourself anymore and you start insulting everything around you"; "It's more like a loss of control. We get a bit carried a way and we just start swearing, whereas in real life we wouldn't"; "Basically, raging is a passing, temporary state, it doesn't last".

Emotional level

Two major emotions were identified in the participants' speeches. Firstly, 79% of the players (n=19, 44 occurrences, i.e. 2.6% of the total occurrences) experience anger, irritation or nervousness during episodes of rage.

"Basically it bubbles up, I mean I can feel it growing, and after a while it feels like I'll explode, I'll really get angry"; "it all depends on the situation, but it really is deep irritation"; "It's being in the heat of excessive anger, and leads to loss of self-control".

Secondly, the presence of frustration is noted in 54% of the players (n=13, 44 occurrences, i.e. 2.6% of total occurrences), for example, "it's frustrating", "I feel very frustrated", "it's going to frustrate me', "it's a cocktail in the end, frustration, irritation, excitement…".

A little under half of the players also describe negative emotions such as guilt or disappointment (n=13, 38 occurrences, i.e. 2% of total occurrences), "It can go really far in terms of emotions when I'm raging, I get to the point where I don't want to do anything, I'm disorientated and annoyed with myself".

Cognitive level

Ninety-two percent of the players (n=22, 84 occurrences, i.e. 5% of total occurrences) express cognitive difficulties. The players describe concentration difficulties: "When I rage, I'm actually less concentrated on the game."; "I have trouble thinking"; "I play without paying much attention"; "I'm not thinking about what I want to do"; "Firstly, I can't concentrate anymore".

The players also describe more impulsivity in their decision-making: "You don't really think about the consequences of these actions straight away"; "well, basically I just don't think! Because I'm spontaneous".

Meanwhile, they describe the presence of external attributions, which manifest in blaming the mistake on others: "we're going to blame it (the mistake) on our teammates"; "we hold them (the other players) a little responsible"; "we're not as kind or lucid"; "I'm moaning at the game in fact".

Twenty one percent of the players (n=5, 13 occurrences, i.e. 0.8% of total occurrences, describe dissociative phenomena, such as: "I'm in another state, in which time goes by faster"; "what we say isn't really relevant to what is actually happening in the game round".

One player describes a vicious circle: "To know we're raging and that we can't control it makes us rage even more than what had happened before, in fact. It's a vicious circle that can send you into a spin, uh… yeah, it can really make us blow up".

Behavioral level

Regarding behavior, all the participants (n=24, 216 occurrences, i.e. 12.9% of total occurrences) mention behavioral manifestations, essentially of aggressiveness. "I speak a little more loudly, and once or twice I've hit my desk but that happens very rarely"; "I'll hit my desk, hit my keyboard, hit the mouse, hit the screen, insult all the players… you know (laughter)"; "I start insulting people and, like, knocking the mouse about, I sit back down, it's true that I can move the mouse abruptly, bang my fist on the table, that kind of thing, but it doesn't go far."; "It can create nervous trembling, a way of acting that highlights aggressiveness in terms of behavior, something that really gives off anger, you can be really agitated, move, and in terms of the tone of voice too".

Verbal aggressiveness is noted for 88% of our participants (n=21, 127 occurrences, i.e. 7.6% of total occurrences). They thus start to "insult other people"; "let go verbally", "I let off steam, I say what I need to say, insults"; "using strong swear words in a spontaneous way".

Physical aggressiveness is present 63% of the players (n=15, 28 occurrences, i.e. 1.7% of total occurrences) and manifests through violence towards pieces of furniture: "Sometimes I might hit or punch the bed"; "throwing my controller"; "hitting the desk"; "for example, throwing objects on the floor. Like knocking a bottle of water that's beside you, and running the risk of damaging your equipment".

Furthermore, 50% of the players (n=12, 63 occurrences, i.e. 4% of total occurrences) can experience a feeling of loss of control: "I'm spontaneous. Whenever I think something, I say it, cash, and I don't think of the consequences it can have";"you just can't control yourself anymore"; "you don't control the situation anymore".

However, 50% of our participants (n=12, 47 occurrences, i.e. 2.8% of total occurrences) add that this phenomenon is a transitory state: "it doesn't last"; "it's not necessarily long" or even "it's temporary".

Causes

All the participants (n=24) evoke external factors to explain the triggering of this phenomenon. Two major causes can be noted in their descriptions. 83% of them (n=20, 133 occurrences, i.e. 7.9% of total occurrences) describe that rage may be caused by interactions with other players: "when others start to insult you'; "when the guy in front of you says 'easy', well, that just drives you crazy"; "Generally, you don't rage when you're alone. You get angry, you'll want to go further in the game, but you won't rage like a guy can make you rage, you know"; "when I play with other people, when the person has a bad attitude, well, I mean when they do badly what they're supposed to do, when they don't follow the logic of the game".

Another cause of rage is experiencing failure, for 88% of the players (n=21, 106 occurrences, i.e. 6.3% of total occurrences), when "we lose many times in a row"; "when we lose"; "the fact of losing"; "a dude that kills me several times in a row"; "when we lose many times in a row"; "when I lose my ranking"; "repeated situations of failure"; "when you just get killed continuously".

Finally, 67% (n=16) of the players link the competitive aspect of video games with the emergence of a rage episode, and describe a greater presence of the phenomenon in "competitive games", when you "measure yourself against the other", and the call for a "competitive spirit"; "even if you're here to play, all of this is still a match so it's a competition"; "people with a very, very competitive spirit mostly. But then I think that everybody is prone to raging if you want to win".

Two older participants (28 and 24 years old) also noted an effect of the player's age.

These gamers explain that in their view, young players are not only a cause of rage, because of their provocative behavior, but are also more likely to experience these episodes of loss of control. "Generally, the kids, those are the ones who are good. They have more reflexes, they're quick, but at the same time, they don't have the maturity to calm down; 'Those who are between 14 and 18 years old, generally those are the ones who fire off comments (aggressive ones) and it's going to piss you off. At the same time, if they lose, they'll blow up faster".

Consequences

For 63% of the players (n=15, 40 occurrences, i.e. 2.4% of total occurrences), rage leads to the loss of the game round. "The round will finish faster and generally, it will end in a defeat"; "my level of gaming worsens"; "it's going to affect the way they play, and they won't play as well and I'm not as likely to win"; "If you're too submerged by rage, yes, you're going to play any old how, and because of that you're going to lose"; "so that makes you lose the round or even the game".

Furthermore, 63% of the players (n=15, 28 occurrences, i.e. 1.7% of total occurrences) explain that rage can have a contagious effect on the other players : "the others also start to get angry and to rage, and then it's over"; "it spreads"; "it'll infect the others"; "it can really destroy the game and its atmosphere"; "when we get into this rage and the others share it and infect the others, and then everybody is angry"; "you can sense the atmosphere, even the players on your team change their way of playing, actually, you can sense it in the movements and the way they play"; "in gaming (online), the problem is that it's basically going to infect the others. It can really spoil the game and its atmosphere, because if everybody is swearing or being sworn at"; "the team's cohesion is challenged, there's no longer any cohesion at all in my view…".

A loss of pleasure in playing or a loss of interest for the game is described by 66% of the players (n=16, 73 occurrences, i.e. 4% of total occurrences). "I don't feel like playing anymore, and I don't enjoy it anymore"; "Either I don"t want to play at all anymore, or I stop playing, or I continue playing but I don't really want to continue"; "It'll also impact my enjoyment, I'll enjoy playing less, you just want to win (…). Actually, you want it to end fast and to move on (laughter)".

Moreover, a positive and liberating effect of these rage episodes is described by 42% of the players (n=10, 36 occurrences, i.e. 1.9% of total occurences). "In the heat of the moment, it clearly did me good. Clearly"; "It's more like an outlet, where everything is allowed basically, and so I'm not going to contain my rage, but offload, for example after a crap week. I'll offload and it'll do me good in fact".

Rage regulation

In order to cope with this phenomenon of rage, the players describe establishing regulation strategies such as changing activities or rationalization.

Changing activities is used by 79% of the players (n=19, 72 occurrences, i.e. 4.3% of total occurrences), for example: "I watch a movie or I read a book"; "I go and get some fresh air"; "Either I do something else, I smoke a cigarette, or I just take a break, I go and eat something". 71% of the players (n=17, 49 occurrences, i.e. 2.9% of total occurrences) try to rationalize their anger. This is conveyed in statements such as: "I manage to tell myself that it's just a game"; "there's no reason to get angry for something like that, there are more important things in life"; "I know how to get some distance on this, I know it's a game and that even if I lose it's not the end of the world"; "Telling yourself that even if you have objectives in the game, just telling yourself that it's only a game round"; "I imagine him (the opponent) as someone inferior to me, to reassure myself in a way".

Lastly, we note that social support may have a role in the regulation of the phenomenon of rage, for 58% of the players (n=14, 40 occurrences, i.e. 2.4% of total occurrences) indicate that "playing with friends, that helps"; "we lose control more easily alone than with friends"; "you rage more when you're alone (in the room)"; "Well, it's the people who are next to me, like my parents, who say 'shut up, stop shouting' (laughter)".

Because personally, I don't realize it".

A majority of the players also described a type of behavior that seems to be closely associated to rage, actually named "rage-quit". "I rag-quit, that means I leave the game.

Then, I turn off the computer directly and I go and relax"; "Sometimes I can leave the game directly"; "Instead of getting angry and insulting the player, I just turned the game off as so not to smash everything"; "what's on repeat in my mind the most is 'do I rage-quit?'. Even if that means getting a warning".

A player also described another strategy which consists in only cutting short the conversations with the others: "I do what we call a mute, basically I turn off the chat".

Discussion

The aim of this study was firstly to explore and define the characteristics of rage in video gaming. The results of the thematic analysis indicate that a specific experience of loss of control exists in the context of online video games. This experience, called rage or tilt by the players, is characterized by a change in the manner of gaming caused by a feeling of frustration during the game. This phenomenon impacts the gamer on emotional, cognitive and behavioral levels.

Rage seems to be brought on by factors external to the game, mainly the interactions with the other players, and in major contests in which the gamer loses. The social and competitive components of these video games may well contribute to experiencing episodes of rage. On the emotional level, rage occurs due to a loss of emotion regulation abilities, resulting in anger and frustration that can involve hatred. After the rage episode, the players described negative feelings such as guilt and self-disappointment. On a cognitive level, the impairment of reasoning abilities results in a difficulty in decision-making and in the implementation of adaptive strategies. We observe an increase in external attribution of the perceived causes of the failure in the game. After the episode, gamers describe a loss of interest and motivation for gaming. Finally, on a behavioral level, gamers highlight significant restlessness and impulsivity, that could be associated with verbal or behavioral aggressive acts. The gamers described two specific characteristics of rage episodes: (1) they are transitory, lasting from a few minutes to a few hours and (2) they are contagious. The rage involves the gamer himself, as well as his teammates and adversaries.

Rage in some aspect could be compared with toxic and aggressive behavior in the context of violent video games. Toxic behavior is defined by anger and frustration generated by a negative event during the game, leading to toxic communication (harmful, contaminated, and disseminated toxic type) [START_REF] Neto | Studying toxic behavior influence and player chat in an online video game[END_REF]. Few studies address the question of toxic behavior [START_REF] Kordyaka | Towards a unified theory of toxic behavior in video games[END_REF]. Compared with the violent nature of the video games, toxic behavior are more related to the type of game (multiplayer), the competitive characteristic of the game and of the gamer [START_REF] Kordyaka | Towards a unified theory of toxic behavior in video games[END_REF][START_REF] Neto | Studying toxic behavior influence and player chat in an online video game[END_REF]. First, rage can lead to toxic and aggressive behavior; especially anger is involved in the relationship between video games and aggressive behavior [START_REF] Morrison | Sex differences in anger states after violent video game play and the relationship of empathy and emotional invalidation as preaggressive predictors in young adults[END_REF][START_REF] Yao | Violent video games exposure and aggression : The role of moral disengagement, anger, hostility, and disinhibition[END_REF][START_REF] Zheng | Online violent video games and online aggressive behavior among Chinese college students : The role of anger rumination and self-control[END_REF].

Then, as toxic behavior, rage is contagious to other players and can impact their psychological well-being [START_REF] Kordyaka | Towards a unified theory of toxic behavior in video games[END_REF][START_REF] Neto | Studying toxic behavior influence and player chat in an online video game[END_REF]. Toxic and aggressive behavior have an influence on the other players and on their psychological state [START_REF] Zsila | Toxic behaviors in online multiplayer games : Prevalence, perception, risk factors of victimization, and psychological consequences[END_REF].

Gamers describe various ways of handling the episodes. First, the most frequently used method, and described as the most efficient, is stopping the game. One of the most specified strategies is "rage-quit", defined as an abrupt, impulsive cessation of the game due to rage. Gamers describe using rationalization in order to set a distance between them and the game. They indicate that social support plays a major role in the regulation of rage episodes, either by limiting its onset or by an incentive from other players to stop the game in order not to pass on their rage to the other members of the team. The most relevant aspect of rage-quit should be the behavior of the player once they quit the game (Wu et al., 2021). Some players calm down, do something else, and others express their emotions sometimes violently.

Emotion regulation is at the core of rage management and is associated with problematic gaming and its consequences on teammates.

Anger and rage imply aggressive behaviors up to and including violence. Gamers are interested in regulating their rage episodes and emotion regulation can be used as a hook for the prevention of the problematic gaming and in e-sport in which rage has to be regulated for performing and being successful. Interestingly, the mediating role of anger rumination between online violent video games and online aggressive behavior can be moderate by selfcontrol [START_REF] Zheng | Online violent video games and online aggressive behavior among Chinese college students : The role of anger rumination and self-control[END_REF] and stress and low performances in the game are a part of a vicious circle leading to toxic behavior [START_REF] Neto | Studying toxic behavior influence and player chat in an online video game[END_REF]. The exposure to repeated toxic behavior during the game, may lead to anxiety, low self-esteem [START_REF] Kordyaka | Towards a unified theory of toxic behavior in video games[END_REF]. As a consequence, helping the player to identify the first signs of a rage episode, informing them that they should stop or change for a less stressful game [START_REF] Neto | Studying toxic behavior influence and player chat in an online video game[END_REF], could be interesting to limit the impact of rage on the player's psychological health, to preserve other players and the game environment [START_REF] Zsila | Toxic behaviors in online multiplayer games : Prevalence, perception, risk factors of victimization, and psychological consequences[END_REF]. Finally, rage could be the visible side of a problematic behavior for the gamer and its relatives (family, game partners). When repetitive and frequent, rage could be a gateway to address problematic gaming and more generally the negative consequences of gaming.

Limits

This study presents several limitations. First, this qualitative research was conducted with a limited sample of 24 online video gamers in three categories. FPS, MOBA, and MMORPG. The population was young and included mostly men, thus excluding a generalization of the results to the general population of online video gamers. Few literature on rage, anger, and aggressive behavior in the context of (violent) video games take into consideration gender. While men and women report an equal level of anger [START_REF] Deffenbacher | State-Trait Anger Theory and the utility of the Trait Anger Scale[END_REF], there is no consensus regarding gender differences in aggression. [START_REF] Fahlgren | Gender Differences in the Relationship between Anger and Aggressive Behavior[END_REF] found that men and women are similarly aggressive but among men aggression is associated with anger trait [START_REF] Fahlgren | Gender Differences in the Relationship between Anger and Aggressive Behavior[END_REF]. In the video game framework, women reported a higher level of anger state than men after playing violent video games [START_REF] Morrison | Sex differences in anger states after violent video game play and the relationship of empathy and emotional invalidation as preaggressive predictors in young adults[END_REF]. Considering the scare literature, more studies is needed to explore gender issues, women may have a different definition, intensity in their rage episodes, and may managed the episode using different strategies from men.

Most participant (66,5 %; n=16) reported a problematic use of video games, which makes it difficult to generalize to the entire population of players. This prevalence is very high, considering we did not specifically look for participants with a problematic use. Rage description may concern more specifically players with a potential problematic use of video game. In addition, for each participant, data from the online questionnaire could not be matched with the data from the interviews (participants did not provide an anonymous code).

It could have been relevant to link the description of the rage and the level of problematic gaming to identify differences in the experience and frequency of rage. This should be included in future studies on rage and on how it impacts gamers behaviors.

In addition, the interviews were conducted via Skype. This approach does not have the advantage of face-to-face interviews, but was adapted to the methodology in order to foster the diversity of the gamers involved in the study. Finally, a qualitative analysis involves a degree of subjectivity, that we tried to limit by using two investigators to evaluate each interview.

Conclusions

This exploratory study is the first to focus on an experience which is specific to online video games among adults. Raging in video games could be defined as a behavioral, emotional and cognitive loss of control, that may spread to other gamers through online interactions. Futures studies should focus on the quantitative assessment of rage and its relations with online gaming disorder, gaming toxic behaviors and its consequences on the psychological well-being of the gamer and his team.

Author Disclosure

Role of Funding Sources

No funding sources

Contributors

AM, EC and HC designed the study and wrote the protocol. AB, VP, MT et JM collected the data and conducted the analysis. AM and EC drafted the first version of the manuscript. H.

Chabrol participated in the study design and drafting the manuscript. All authors edited and approved the final manuscript.

Conflict of interest

Axelle Moreau, Amelie Bethencourt, Valentin Payet, Mégane Turina, Julia Moulinard, Henri Chabrol, and Emeline Chauchard declare no conflict of interest.