Monireh Dehkhodaei

Mehdi Sahihi

Amiri Hadi

Sajjad Rudbari

Reza Gharaghani

Salman Azadbakht

Abolghasem Abbasi Taheri

Kajani

Hadi Amiri Rudbari

Sajjad Gharaghani

Reza Azadbakht

Salman Taheri

Abbasi Kajani

Studies of DNA-and HSA-binding properties of new nano-scale green synthesized Ni (II) complex as anticancer agent using spectroscopic methods, viscosity measurement, molecular docking, MD simulation and QM/MM

Keywords:

Introduction

More than 40% of compounds with pharmacological activity that are introduced as drugs have poor water solubility problem. Beside, reducing the size of drugs to nano-scale causes the several unique properties that could not be seen in their bulk form; including high surface to volume ratio, unique magnetic behavior and adjusting of their dispersibility in various solvents [START_REF] Gautam | PAMAM dendrimers: novel polymeric nanoarchitectures for solubility enhancement of candesartan cilexetil[END_REF]. Hence, synthesis of nano-scale drugs has attracted many researchers due to the abovementioned properties [START_REF] Dreaden | The golden age: gold nanoparticles for biomedicine[END_REF], [START_REF] Jain | Gold nanoparticles as novel agents for cancer therapy[END_REF], [START_REF] Parveen | Nanoparticles: a boon to drug delivery, therapeutics, diagnostics and imaging, Nanomed[END_REF]. Schiff base compounds are considered as a very important class of organic materials, due to their interesting chemical and physical properties [START_REF] Panyam | Biodegradable nanoparticles for drug and gene delivery to cells and tissue[END_REF]. The Schiff base ligands have high coordination capability and many strong biological properties. Although these compounds have been extensively studied, they are still highly regarded [START_REF] Keypour | Synthesis of two new N2O4 macroacyclic Schiff base ligands and their mononuclear complexes: spectral, X-ray crystal structural, antibacterial and DNA cleavage activity[END_REF]. Schiff base ligands exhibited the biological activity such as antitumor, antibiotics, antifungal, antimicrobial, and anticancer properties [START_REF] Abo-Aly | Spectroscopic and structural studies of the Schiff base 3-methoxy-N-salicylidene-o-amino phenol complexes with some transition metal ions and their antibacterial, antifungal activities[END_REF], [START_REF] El-Wakiel | Synthesis, spectral, antitumor, antioxidant and antimicrobial studies on Cu (II), Ni (II) and Co (II) complexes of 4-[(1Hbenzoimidazol-2-ylimino)-methyl]-benzene-1, 3-diol[END_REF], [START_REF] Mahmoud | Novel Schiff base ligand and itsmetal complexes with some transition elements. Synthesis, spectroscopic, thermal analysis, antimicrobial and in vitro anticancer activity[END_REF], [START_REF] Rahman | Nano structure iron (II) and copper (II) Schiff base complexes of a NNO-tridentate ligand as new antibiotic agents: spectral, thermal behaviors and DNA binding ability[END_REF]. Also, they have appropriate affinity for binding to HSA and DNA [START_REF] Kazemi | Synthesis, characterization and separation of chiral and achiral diastereomers of Schiff base Pd (II) complex: a comparative study of their DNA-and HSA-binding[END_REF], [START_REF] Kazemi | Synthesis, characterization and biological application of four novel metal-Schiff base complexes derived from allylamine and their interactions with human serum albumin: experimental, molecular docking and ONIOM computational study[END_REF], [START_REF] Khosravi | Synthesis, characterization, crystal structure and HSA binding of two new N, O, O-donor Schiff-base ligands derived from dihydroxybenzaldehyde and tert-butylamine[END_REF]. The azomethine linkage in Schiff bases is responsible for their biological activity and several studies proved that the presence of lone pair electrons in sp 2 hybridized orbital of nitrogen atom of the azomethine group plays a remarkable chemical and biological role. On the other hand, the metal ions present in complexes, not only accelerate the drug activity of compounds, but also increase the effectiveness and efficiency of the ligands [START_REF] Anacona | Synthesis, characterization and antibacterial activity of a Schiff base derived from cephalexin and sulphathiazole and its transition metal complexes[END_REF], [START_REF] Mohini | Synthesis, antimicrobial and anti-biofilm activities of novel Schiff base analogues derived from methyl-12-aminooctadec-9-enoate[END_REF]. Hence, great effort has been devoted to design and development of Schiff base complexes using different metal ions, which possess excellent antioxidant, antibacterial and antitumor activities. Nickel plays a role in urease as a member of redox enzymes family [START_REF] Eitinger | Nickel transport systems in microorganisms[END_REF], [START_REF] Ragsdale | Nickel biochemistry[END_REF]. Also, considerable biological activities are reported for Ni(II) Schiff base complexes [START_REF] Kazemi | Synthesis, characterization and biological application of four novel metal-Schiff base complexes derived from allylamine and their interactions with human serum albumin: experimental, molecular docking and ONIOM computational study[END_REF], [START_REF] Zangrando | Synthesis, characterization and bio-activity of nickel (II) and copper (II) complexes of a bidentate NS Schiff base of S-benzyl dithiocarbazate[END_REF]. It is clear that many compounds utilize their drug effects through binding to DNA or carrier proteins such as Human Serum Albumin (HSA). Studying on the interactions of these metal complexes with DNA or HSA is critical to design target specific, more efficient and less toxic drugs. Interaction of candidate drugs with these biomacromolecules could be reflected using spectroscopic (fluorescence, UV-Vis and so on) and computational methods (molecular docking, molecular dynamics simulation and qm/mm). Recently, we have focused on designing and synthesis of novel Schiff base ligands, and investigation of their assembly with metals as well as their biological activities [START_REF] Kazemi | Synthesis, characterization and separation of chiral and achiral diastereomers of Schiff base Pd (II) complex: a comparative study of their DNA-and HSA-binding[END_REF], [START_REF] Kazemi | Synthesis, characterization and biological application of four novel metal-Schiff base complexes derived from allylamine and their interactions with human serum albumin: experimental, molecular docking and ONIOM computational study[END_REF], [START_REF] Khosravi | Synthesis, characterization, crystal structure and HSA binding of two new N, O, O-donor Schiff-base ligands derived from dihydroxybenzaldehyde and tert-butylamine[END_REF]. Herein, a new Schiff base Ni(II) complex (NiL2) was synthesized in water as a green solvent and at ambient temperature. These conditions are valuable to introduce a compound as a drug candidate. Then the compound was characterized using FT-IR and elemental analysis. In addition, its molecular structure was determined by single crystal X-ray diffraction technique. The cell viability percent of HeLa cancer cells was first studied by MTT assay. In order to increase the colloidal stability and suitability for the biomedical applications, the nano-scale compound was also synthesized and used for in vitro studies. Finally, binding ability of the nano-scale Ni(II) Schiff base complex with calf thymus DNA (CT-DNA) and HSA was investigated using combination of experimental (UV-Vis, fluorescence, circular dichroism (CD) and viscosity) and computational (molecular docking, molecular dynamics simulation and qm/mm) methods.

Experimental section

Materials sources

All of the used chemicals for synthesis of the complex including 3-amino-prop-1-ene (Allylamine), 2-hydroxybenzaldehyde (salicylaldehyde), triethylamine and nickel acetate (Ni(OAc)2) were purchased from Merck Co. and were used without further purification. All of the used salts for buffer preparation were analytical grade and were dissolved in double distilled water. All of the solutions were used freshly after preparation. Also, human serum albumin (HSA), CT-DNA, RPMI-1640 medium, Fetal bovine serum (FBS), dimethyl sulfoxide (DMSO), antibiotics (penicillin-streptomycin) solution, and 3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide (MTT) were obtained from Sigma-Aldrich.

Synthesis of the bulk and nano-scale Ni (II) Schiff base complex (NiL2)

The solution (5 ml) of Allylamine (IUPAC name: 3-amino-prop-1-ene) (4 mmol) in double distilled water was added slowly to 5 ml of a methanolic stirred solution of salicylaldehyde (IUPAC name: 2-hydroxybenzaldehyde) (4 mmol) in ambient temperature. The color immediately changed to yellow and the mixture was then stirred for 2 h before removal of solvent under vacuum. A solution of triethylamine (4 ml) in double distilled water was added dropwise to ligand solution. The mixture was stirred for 15 min again. Then, a solution of appropriate metal salt (Ni (OAc)2 (2 mmol)) in absolute deionized water (60 ml) was added to the mixture, gently, for synthesis of the complex in bulk. The resulting solution was stirred for 1 h in room temperature. After concentration, dark green precipitate was collected by filtration and was washed with methanol and deionized water several times. Appropriate single crystals for Xray crystallography were obtained directly from the reaction mixture. In order to prepare the complex in nano-scale; ultrasonic probe was fixed in ligand solution and was positioned under the high-intensity ultrasonic waves. At this time, a solution of Ni (OAc)2 (2 mmol) in 60 ml deionized water was added to the mixture gently under ultrasonic condition for 2 h. The green powder was isolated from the solution and was purified by washing with methanol and deionized water several times and was allowed to get dried at room temperature.

Single crystal diffraction studies

X-ray data for NiL2 were collected on a STOE IPDS-II diffractometer with graphite monochromated Mo-Kα radiation. For NiL2, a green crystal was chosen using a polarizing microscope and was mounted on a glass fiber which was used for data collection. Data were collected at temperature of 298 K in a series of [START_REF] Gautam | PAMAM dendrimers: novel polymeric nanoarchitectures for solubility enhancement of candesartan cilexetil[END_REF] scans in 1° oscillations and were integrated using the Stöe X-AREA [START_REF] Stoe | X-AREA, Version 1.30, Programfor the Acquisition and Analysis of Data[END_REF] software package. A numerical absorption correction was applied using the X-RED [START_REF] Stoe | X-RED: Programfor Data Reduction and Absorption Correction[END_REF] and X-SHAPE [START_REF] Stoe | X-SHAPE, Version 2.05, Program for Crystal Optimization for Numerical Absorption Correction[END_REF] softwares. The data were corrected for Lorentz and Polarizing effects. The structure was solved by direct methods using SIR2004 [START_REF] Burla | SIR2004: an improved tool for crystal structure determination and refinement[END_REF]. The non-hydrogen atoms were refined anisotropically by the full-matrix least-squares method on F 2 using SHELXL [START_REF] Sheldrick | Crystal structure refinement with SHELXL[END_REF]. All hydrogen atoms were added at ideal positions and were constrained to ride on their parent atoms.

Cell viability assay

The potential anticancer activities of the bulk solution and nano-scale NiL2 dispersion were studied using MTT assay on human breast cancer cell line ((MCF-7)) according to the previously reported procedure [START_REF] Kajani | Gold nanoparticles as potent anticancer agent: green synthesis, characterization, and in vitro study[END_REF]. Briefly, the cancer cells were first cultured in RPMI-1640 medium supplemented with 10% FBS and 1% antibiotics solution and were maintained in a humidified 5% CO2 incubator at 37 °C for two weeks. After propagation of the cells and obtaining the appropriate density, the cells were harvested and were seeded on 96-well plates at a density of 104 cells per well containing 200 μl medium. The cells were subsequently incubated overnight at the same conditions before treatment with different concentrations [START_REF] Rahman | Nano structure iron (II) and copper (II) Schiff base complexes of a NNO-tridentate ligand as new antibiotic agents: spectral, thermal behaviors and DNA binding ability[END_REF][START_REF] Stoe | X-RED: Programfor Data Reduction and Absorption Correction[END_REF][START_REF] Lakowicz | Plasmonics in biology and plasmon-controlled fluorescence[END_REF] and 100 μM) of the bulk and nano-scale NiL2 complex for 48 h. The medium was then removed and 100 μl MTT solution (0.5 mg/ml in media) was added into each well and the plates were incubated again at 37 °C for 4 h. Finally, the medium was carefully discarded and the remained formazan crystals were dissolved in 150 μl DMSO and the absorbance was measured at 570 nm. Three independent experiments were conducted for each toxicity endpoint and the results were presented as the mean values obtained from three independent experiments. The cell viability was determined as ratio of absorbance values from each treatment and the control. To further evaluate the anticancer effects of the compounds, the cell morphology changes after exposure to the compounds were also monitored by optical microscope.

Preparation of DNA, HSA and nano-scale NiL2 for binding experiments

The stock solution of CT-DNA was prepared in 50 mM Tris buffer at pH 7.5 using double-distilled deionized water and was stored at 4 °C. The CT-DNA concentration per nucleotide was determined using absorption intensity at 260 nm after adequate dilution with the buffer and using the reported molar absorptivity of 6600 M -1 •cm -1 [START_REF] Reichmann | A further examination of the molecular weight and size of desoxypentose nucleic acid[END_REF]. Purity of CT-DNA solution was confirmed by ratio of UV absorbance at 260 and 280 nm (A260/A280 = 1.9), indicating that CT-DNA is free from protein impurity [START_REF] Zheng | Synthesis, structure and molecular docking studies of dicopper (II) complexes bridged by N-phenolato-N′-[2-(dimethylamino) ethyl] oxamide: the influence of terminal ligands on cytotoxicity and reactivity towards DNA and protein BSA[END_REF]. Also, a stock solution of HSA was prepared by dissolving the desired amount of HSA in 50 mM phosphate buffer (pH = 7). The HSA stock solution was stored at 4 °C in the dark and was used within 2 h. HSA concentration was determined by UV-Vis spectrophotometry using the molar absorption coefficient 35,700 M -1 •cm -1 at 278 nm [START_REF] Pace | How to measure and predict the molar absorption coefficient of a protein[END_REF]. The appropriate amount of the nano-scale NiL2 was dispersed in Tris and phosphate buffer for DNA and HSA binding experiments, respectively.

Fluorescence spectroscopy measurements

Fluorescence quenching experiments were carried out using quartz cuvette with 1 cm optical path length and the excitation and emission slits were set at 5 and 10 nm, respectively. In our primary experiments, ethidium bromide (EB) emission was checked in the presence of various amounts of CT-DNA. The results showed that the emission of EB was increased up to mole ratio of CT-DNA:EB = 10:1 and there was no significant increasing in emission after the mentioned mole ratio. Hence, the CT-DNA solution was stirred with EB with molar ratio of CT-DNA:EB 10:1 and was incubated for 1 h at 4 °C for completion of interaction between CT-DNA and EB. Then, various amounts of nano-scale NiL2 dispersion (250 μM) were added to the mixture of CT-DNA:EB. The fluorescence spectra were measured in the range of 500-700 nm with exciting wavelength at 520 nm. In each measurement after addition of NiL2, the mixture was allowed to stand for 2 min. Moreover, the measured fluorescence intensities were corrected for the dilution and the inner-filter effect. To eliminate the inner filter effects, absorption measurements were carried out at the fluorescence excitation and emission wavelengths. The extent of this effect can be roughly evaluated with the following relationship: [START_REF] Bi | Binding of several anti-tumor drugs to bovine serum albumin: fluorescence study[END_REF] where, Fcorr and Fobs are the corrected and observed fluorescence intensities, respectively, and Aex and Aem are the absorption of NiL2 at excitation and emission wavelengths, respectively. All the fluorescence intensities obtained in this study were corrected. For investigation of the HSA interaction with NiL2 by this technique 2 ml of HSA solution (10 μM) was placed into the quartz cell and various amounts of NiL2 (200 μM) were added to the cell. The fluorescence intensities were measured with excitation wavelength at 295 nm and emission wavelength range of 300-450 nm. In each measurement, the mixture was allowed to incubate for 2 min after addition of the complex. Furthermore, all intensities were corrected for the dilution and inner filter effect in fluorescence experiments.

UV-Vis absorption spectroscopy measurements

Absorption spectroscopy titration experiment was performed by addition of various amounts of CT-DNA (1 mM) to 2 ml of nano-scale NiL2 dispersion (200 μM). All NiL2-CT-DNA solutions were allowed to incubate for 2 min before recording the related spectra. Absorption curves of NiL2-CT-DNA mixtures were corrected by subtracting the spectra of CT-DNA and all intensities were corrected for the dilution effect. The UV-Vis absorption spectra of HSA solution (10 μM) in the absence and presence of various amounts of the nano-scale NiL2 dispersion (200 μM) were recorded after 2 min of incubation. All intensities were corrected for the dilution effect.

Circular dichroism spectroscopy

The circular dichroism (CD) spectra were recorded at room temperature using 0.1 cm quartz cell in the far-UV region (200-260 nm) for HSA and in the far and near-UV (205-310) for DNA to clarify the structural changes of biomacromolecules during their interaction with NiL2. The CD spectrum of HSA and DNA solution was recorded before and after addition of the nano-scale NiL2 complex with molar ratio of 1:1.

Viscosity measurement

Viscosity experiments were carried out using a rotational viscometer and the measurements were performed at 200 rpm at room temperature. The viscosity of CT-DNA solution was measured in the presence of various amounts of the nano-scale NiL2 dispersion. The obtained data are presented as (η/η 0 ) 1/3 versus [NiL2] / [CT-DNA], where η and η 0 are the viscosity of CT-DNA in the absence and presence of NiL2, respectively.

Computational section

Molecular docking calculations

Docking study was carried out to indicate the HSA and DNA-binding site for NiL2 complex. The 3D structure of NiL2 was obtained using the .cif file of its X-ray crystal structure. The .cif file was converted to the .pdb format using the Mercury software (http://www.ccdc.cam.ac.uk/). The crystal structures of HSA (PDB ID: 1AO6) and DNA (PDB ID: 423D) with sequence d(ACCGACGTCGGT)2 were taken from the Brookhaven Protein Data Bank (http://www.rcsb.org/pdb). The resolution of these files was 2. 5 and 1.6 Å for HSA and DNA, respectively. Water molecules of the .pdb files were deleted, missing hydrogen atoms and Gasteiger charges were added. Flexible ligand docking was carried out by AutoDock 4.2.5.1 molecular docking program using the implemented empirical free energy function and the Lamarckian Genetic Algorithm [START_REF] Morris | Automated docking using a Lamarckian genetic algorithm and an empirical binding free energy function[END_REF]. The Gasteiger charges were added to prepare the macromolecule input file for docking and the AutoGrid was used to calculate grids. For docking of the NiL2 complex to HSA the grid box was centred on Cα of the Trp-214 residue of protein and a docking with 90 lattice points along X, Y, and Z axes was performed to find the active site of complex to the HSA. Grid point spacing was set to 0.375 Å, to allow the complex to rotate freely. For the docking of NiL2 complex with DNA, a blind docking with 126 lattice points along X, Y, and Z axes was performed to find the binding site of complex on DNA with a grid point spacing of 0.375 Å, to allow the complex to rotate freely. In the next step, the centre of the grid box was located at the binding site and the second docking was performed using a cubic box with 90 × 90 × 90 dimensions. 250 docking runs with 25,000,000 energy evaluations for all runs were performed.

Molecular dynamics (MD) simulation

MD simulation method was used to compare the structural changes of DNA and HSA in the absence and presence of NiL2 complex. The NiL2-HSA and NiL2-DNA complexes with the most negative free binding energy were considered as the initial conformations for the MD studies. GROMACS 5.1.2 package and Amber99 force field [START_REF] Pérez | Refinement of the AMBER force field for nucleic acids: improving the description of α/γ conformers[END_REF], [START_REF] Wang | Development and testing of a general amber force field[END_REF] were used to carry out all MD studies. The partial atomic charges of NiL2 were calculated using Gaussian 03 at the level of B3LYP/6-31G** [START_REF] Lee | Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density[END_REF], while other intramolecular force-field parameters were generated with ACPYPE [START_REF] Da Silva | ACPYPE-antechamber python parser interface[END_REF] and VFFDT [START_REF] Zheng | VFFDT: a new software for preparing AMBER force field parameters for metal-containing molecular systems[END_REF], respectively for the organic moiety and the metal center. Free biomacromolecules and NiL2-biomacromolecule complexes were located in the cubic box with the periodic boundary conditions in the three directions. The solutes were placed in the center of box and the minimum distance between solute surface and the box was 1.0 nm. The box filled with TIP3P water molecules [START_REF] Mark | Structure and dynamics of the TIP3P, SPC, and SPC/E water models at 298 K[END_REF], [START_REF] Jorgensen | Comparison of simple potential functions for simulating liquid water[END_REF], and the solvated systems were neutralized by adding appropriate amounts of sodium ions (Na + ) and chloride ions (Cl -). After energy minimization using the steepest descent method, the systems were equilibrated for 200 ps at the temperature of 300 K. Finally, a 10 ns MD simulation was carried out at 1 bar and 300 K. Parrinello-Rahman barostate [START_REF] Parrinello | Polymorphic transitions in single crystals: a new molecular dynamics method[END_REF], [START_REF] Rahman | Molecular dynamics study of liquid water[END_REF] at 1 bar, Berendsen thermostat [START_REF] Lemak | On the Berendsen thermostat[END_REF] at 300 K, 9 Å cut off for van der Walls and Coulomb interactions and the particle mesh Ewald (PME) method [START_REF] Darden | Particlemesh Ewald: an N•log (N)method for Ewald sums in large systems[END_REF], [START_REF] Essmann | A smooth particle mesh Ewald method[END_REF] for long range electrostatics were used. The leap-frog algorithm with the 2 fs time step was used to integrate the equation of motions. Finally, an all-bond constrain was used to keep the NiL2 from drifting and the atomic coordinates were recorded to the trajectory file every 0.5 ps for later analysis.

QM/MM study

Own N-layered Integrated Molecular Orbital and molecular Mechanics (ONIOM) methodology was employed to carry out QM/MM calculations. The ONIOM can be considered as a hybrid method, quantum mechanical method (QM) with a molecular mechanics (MM). This method is able to apply different ab initio or semi-empirical methods to different parts of a system and combine to produce reliable geometry and energy at reduced computational time [START_REF] Farrokhpour | A rare dihydroxo copper (II) complex with ciprofloxacin; a combined experimental and ONIOM computational study of the interaction of the complex with DNA and BSA[END_REF]. In the two layers ONIOM method, the total energy (EONIOM) of the whole system is a sum of the three independent energy calculations: [START_REF] Svensson | ONIOM: amultilayered integrated MO+MMmethod for geometry optimizations and single point energy predictions. A test for Diels-Alder reactions and Pt (P (t-Bu) 3) 2+ H2 oxidative addition[END_REF], [START_REF] Vreven | Combining quantum mechanics methods with molecular mechanics methods in ONIOM[END_REF] Real system contains full geometry of the molecule and is considered as MM layer while the model system contains the chemically most important (core) part of the system that is considered as QM layer. Furthermore, the stability of respective adducts was obtained by the evaluation the interaction energy (ΔE), which is calculated by the following equation:

Where ΔE Biomacromolecule/NiL2, ΔEBiomacromolecule and ΔENiL2 are the energy of the optimized adduct of Biomacromolecule-NiL2, Biomacromolecule and NiL2, respectively. Biomacromolecule is DNA and HSA here as receptor. A two layers QM/MM method was used for all calculations. HSA and DNA are considered as low layer while the NiL2 is considered as high layer in the current work. Molecular mechanics method (UFF) and semi empirical method (PM6) were selected for the low layer and the high layer, respectively.

Results

Analysis of synthesis and single crystal X-ray diffraction study

The reaction of Schiff base ligand with metal ion is shown in Scheme 1A. This complex was synthesized in water as solvent that is an important advantage for a drug candidate. The complex was obtained in good yield (more than 80%). The structure of complex was confirmed by Singlecrystal X-ray diffraction and FT-IR techniques. The result of X-ray crystallography (ORTEP representation) is shown in Scheme 1B. The crystallographic data and selected bond lengths and angles are collected in Tables S1 andS2. In the solid state, NiL2 is placed in a crystallographic center of symmetry that is right in the middle point of the Ni(II) in a planar-transoid conformation (Scheme 1B). The molecular unit is centrosymmetric and is made of equivalent halves. The crystallographic data reveal that the NiL2 complex is four-coordinated by two phenolate oxygen and two imine nitrogen atoms of two Schiff base ligands. The Schiff base ligands are arranged with respect to each other in a trans orientation such that one allyl substituent lies above the Although there are no classical hydrogen bonds in the structure, but two intramolecular C H ⋯ O contacts are observed in asymmetric unit of structure (Scheme 1B). These no classical hydrogen bonds and also distances and angles of these interactions are similar to those are seen in related complexes [START_REF] Kazemi | Synthesis, characterization and separation of chiral and achiral diastereomers of Schiff base Pd (II) complex: a comparative study of their DNA-and HSA-binding[END_REF]. Also, FT-IR spectrum of NiL2 complex is shown in Fig. S1.

Dynamic light scattering (DLS) and scanning electron microscopy (SEM) measurements

Fig. 1A shows the size distributions of NiL2 dispersion. Three groups, one corresponding to very small particles (2-7 nm), one corresponding to small particles (100-300 nm) and the other representing large aggregates (5000-6000 nm), are observed in DLS distribution diagram. As it can be seen in this figure, about the 96% of synthesized particles have the size under 300 nm. Also, the morphology of synthesized nanoparticles was investigated under SEM. The particles exhibited an approximately plate shape (Fig. 1B).

MTT assay

The anticancer activity was considered as one of the most important potential applications of the compounds. To this aim, the cell viability percent of HeLa cancer cells after 48 h exposure to different concentrations of NiL2 solution was first studied by MTT assay. The results (Fig. 2) are indicator of the dose dependent anticancer activity, as the cell mortality increased in the higher concentrations of the compound. In the best condition, more than 50% (53.59%) cell mortality was observed after 48 h exposure to 100 μM of the NiL2 solution. Regarding to the obtained hopeful results for the NiL2 solution and also, in order to increase the colloidal stability and suitability for the biomedical applications, the nano-scale NiL2 was also synthesized and used for in vitro studies. The nano-scale compound could also represent the higher anticancer effects due to their probable increased bioavailability. About 20% differences for anticancer activity were observed between the NiL2 solution and nano-scale compound dispersion. The maximum mortality of 79.36% HeLa cells obtained after 48 h exposure to 100 μM of nano-scale compound dispersion (Fig. 2). This result indicated that the anticancer activity of the compound is affected by its size. Regarding to the appropriate anticancer activity of the compound, further studies could be suggested for development of the potent anticancer agents.

Fluorescence study for interaction of nano-scale NiL2 dispersion with CT-DNA and HSA

The fluorescence measurement as sensitive and effective technique for investigation the binding ability of nano-scale NiL2 dispersion with CT-DNA and HSA was performed. Fig. 3A represents the fluorescence quenching of DNA-Ethidium Bromide (EB) spectrum by adding nano-scale NiL2 dispersion. Hydrogen transfer from one of the amino groups of EB into the solution causes non-radiative decay and a weak fluorescence emission. Although, a significant increase in fluorescence intensity of EB is observed in the presence of DNA due to intercalation of the EB molecules into the double helix of DNA [START_REF] Fromherz | Photoinduced electron transfer in DNA matrix from intercalated ethidium to condensed methylviologen[END_REF], [START_REF] Waring | Complex formation between ethidium bromide and nucleic acids[END_REF]. NiL2 can displace EB in a non-competitive manner by changing the DNA conformation. The DNA-bound EB molecules are converted to their free form in solution and causes fluorescence quenching. On the other hand, the intrinsic fluorescence of HSA is mainly due to the tryptophan residue [START_REF] Ganeshpandian | New ruthenium (II) arene complexes of anthracenyl-appended diazacycloalkanes: effect of ligand intercalation and hydrophobicity on DNA and protein binding and cleavage and cytotoxicity[END_REF]. The fluorescence intensity of protein was quenched through the addition of the nano-scale NiL2 dispersion. This implies that the compound strongly interacts with HSA, leading to changes of microenvironment around the Trp-214 residue in HSA [START_REF] Mrkalić | Synthesis of novel palladium (II) complexes with oxalic acid diamide derivatives and their interaction with nucleosides and proteins. Structural, solution, and computational study[END_REF], [START_REF] Wang | Interaction of the flavonoid hesperidin with bovine serum albumin: a fluorescence quenching study[END_REF]. Fig. 4A represents the fluorescence quenching of HSA at the presence of various amount of the compound. The Scachard equation was used to obtain the binding affinity of nano-scale NiL2 dispersion to DNA and HSA: [START_REF] Lakowicz | Plasmonics in biology and plasmon-controlled fluorescence[END_REF] where, F0 and F are the fluorescence intensity of biomacromolecule in the absence and presence of the compound, respectively. [Q] Is the concentration of quencher that quencher is nano-scale NiL2 dispersion here. "Kb" is binding constant and is obtained from the plot of Ln((F0 -F)/F) versus Ln [Q] as y-intercept. These plots are shown in Figs. 3B and4B for DNA and HSA respectively. "n" is the number of binding site per biomacromolecule and is determined as slope of the plot. The values of n are nearly 1 for NiL2 with both of biomacromolecules, indicating that the compound binds to DNA and HSA with molar ratio of 1:1. The calculated results are shown in Table 1. Table 1. The DNA and HSA-binding constants (Kb), the number of binding sites (n), the Stern-Volmer constants (Ksv) and the quenching rate constants (kq) of the NiL2.

Type of biomolecule Kb(M -1 ) (fluorescence) Kb(M -1 ) (UV-Vis) n

KSV(M -1 ) kq(M -1 To determine the interaction mode of the nano-scale NiL2 dispersion with CT-DNA and HSA, the Stern-Volmer quenching plot was obtained by monitoring the fluorescence quenching of biomacromolecules with increasing the concentration of the compound according to the Stern-Volmer equation: [START_REF] Dighe | Synthesis of β-carboline-based n-heterocyclic carbenes and their antiproliferative and antimetastatic activities against human breast cancer cells[END_REF] where, F0 and F are the fluorescence intensity of biomacromolecules in the absence and presence of the compound, respectively. Ksv is the Stern-Volmer quenching constant, kq is the quenching rate constant of biomacromolecules and τ is the average lifetime of biomacromolecules without quencher which is typically equal to 10 -8 s for biomacromolecules. Ksv is determined from the plot of F0/F vs. [Q] (Figs. S2 andS3). Also, the values of Ksv were obtained for both of biomacromolecules and are presented in Table 1. Fluorescence quenching is classified as two mechanisms: static quenching and dynamic quenching. In the static mechanism, the fluorophore and the quencher collide together in the ground state while fluorophore and quencher collide together in the excited state in dynamic mechanism. Linearity of the Stern-Volmer plots indicates that fluorescence quenching has only one mechanism: dynamic or static [START_REF] Lakowicz | Plasmonics in biology and plasmon-controlled fluorescence[END_REF]. In this study, the values of kq were obtained about 10 12 M -1 •S -1 for DNA and HSA. These values are greater than limiting diffusion rate constant of the diffusional quenching for biopolymers (2 × 10 10 M -1 •S -1 ) and support that quenching fluorescence of DNA and HSA occurs by static mechanism.

Thermodynamic parameters and determination of nature of interaction forces

According to the Ross and Subramanian report, there are four non-covalent effective forces in interaction between compounds and biomacromolecules including hydrophobic, hydrogen, electrostatic and Vander-waals interactions [START_REF] Ross | Thermodynamics of protein association reactions: forces contributing to stability[END_REF]. For determination the type of interactions NiL2 with biomacromolecules, the thermodynamic parameters were calculated using the following equations:

where, R is the gas public constant, T is the experimental temperature, K is binding constant at the corresponding temperature, ΔG°, ΔH° and ΔS° are the free energy changes, the enthalpy changes and the entropy changes, respectively. ΔH° and ΔS° were estimated from the plot LnK vs. 1/T as slop and y-intercept, respectively. Also, ΔG° was obtained using Eq. ( 7). The estimated Values are summarized in Table 2. The positive ΔH° and ΔS° values for both of biomacromolecules indicated that hydrophobic interactions have the main role in binding of nano-scale NiL2 dispersion to CT-DNA and HSA. Also, the negative ΔG° values in all experimental temperatures illustrated that the interactions occurred spontaneously. Plots of LnK vs. 1/T for both of biomacromolecules are present in Figs. S4 andS5. Table 2. The thermodynamic parameters for binding of NiL2 to HSA and CT-DNA, including binding constants at the corresponding temperatures (Kb), free energy change (ΔG°), enthalpy change (ΔH°) and entropy change (ΔS°).

Type of biomolecule

UV-Vis absorption

Absorption titration experiment as an operational and very easy method was also carried out for more investigation of HSA and DNA binding of NiL2. Generally, hyperchromic or hypochromic effect and red or blue shift are observed in the UV-Vis spectrum of a drug upon its DNA-binding. Hypochromic and red shift are indicative of intercalation mode involving an interaction between π*-orbital of drug with π-orbital of DNA base pairs [START_REF] Sirajuddin | Drug-DNA interactions and their study by UV-visible, fluorescence spectroscopies and cyclic voltametry[END_REF]. Therefore, energy level of π*-orbital of drug decreases which causes red shift in its UV-Vis spectrum. Furthermore, the coupled π*-orbital is filled and so the probability of electron transition is decreased and hypochromic is observed [START_REF] Sirajuddin | Drug-DNA interactions and their study by UV-visible, fluorescence spectroscopies and cyclic voltametry[END_REF], [START_REF] Abdi | Mono-and polynuclear copper (II) complexes based on 2, 4, 6-tris (2-pyridyl)-1, 3, 5-triazine and its hydrolyzed form[END_REF]. While a groove binding or electrostatic interaction leads to hyperchromic effect along with blue shift [START_REF] Kathiresan | Mixed-ligand copper (II)-phenolate complexes: structure and studies on DNA/protein binding profiles, DNA cleavage, molecular docking and cytotoxicity[END_REF].

The absorption titration was carried out by adding various amount of CT-DNA to nano-scale NiL2 dispersion. The absorption spectrum of nano-scale NiL2 dispersion indicates two absorption bands at 250 and 320 nm that attributed to π → π* and n → π* transitions, respectively (Fig. 3C).

The band in 320 nm was chosen to investigate the interaction of nano-scale NiL2 dispersion with DNA. Although, the existence of hypochromic effect in UV-Vis spectrum attributed to intercalation mode, some of reports predicate that hyperchromic effect is relevant to groove binding due to the changes in solvent distribution and the orientation effect in moving into the more hydrophobic environment of groove [START_REF] Kathiresan | Mixed-ligand copper (II)-phenolate complexes: structure and studies on DNA/protein binding profiles, DNA cleavage, molecular docking and cytotoxicity[END_REF]. In order to assess the binding ability of the compound to CT-DNA, the binding constant (Kb) was estimated by monitoring the changes of absorbance with increasing concentration of CT-DNA and using the following equation: [START_REF] Ozluer | In vitro DNA binding studies of anticancer drug idarubicin using spectroscopic techniques[END_REF] where, εf and εa are the molar absorbance coefficient of NiL2 in the absence and presence of CT-DNA, respectively. εb is the obtained molar absorbance coefficient in saturation and [DNA] is the concentration of CT-DNA. The plot of 1/(εa -εf) versus 1/[DNA] gives Kb as ratio of y-intercept to slop (Fig. 3D and Table 1). Also, the absorption titration was carried out to investigate the interaction of HSA and nano-scale NiL2 dispersion. Fig. 4C shows UV-Vis spectra of HSA in the absence and presence of nano-scale NiL2 dispersion. There is one peak at 278 nm is assigned to π → π* transition of the phenyl rings in aromatic amino acids (Trp, Tyr and Phe) [START_REF] Kathiresan | Mixed-ligand copper (II)-phenolate complexes: structure and studies on DNA/protein binding profiles, DNA cleavage, molecular docking and cytotoxicity[END_REF]. By addition of the nano-scale NiL2 dispersion, a hypochromic effect is observed in HSA spectrum. This result shows that the aromatic amino acids of HSA were exposed to an aqueous environment upon binding to the compound. Following equation was used to investigate the binding affinity of NiL2 to HSA at 278 nm: [START_REF] Lown | Molecular mechanism of binding of pyrrolo (1, 4) benzodiazepine antitumour agents to deoxyribonucleic acid: anthramycin and tomaymycin[END_REF] where, A0 and A are the absorbance of HSA in the absence and presence of the compound, respectively. Amax is the obtained absorbance at saturation and [M] is the concentration of nanoscale NiL2 dispersion. The plot of 1 / (A -A0) versus 1 / [M] gives Kb as ratio of y-intercept to slope (Fig. 4D and Table 1). Generally, binding constant of a drug with a carrier protein such as HSA should be high enough to bind and transfer through body. Moreover, in order to release a drug since arrival at its target, binding constant should not be too high [START_REF] Kyropoulou | Ni (II) complexes with non-steroidal anti-inflammatory drug diclofenac: structure and interaction with DNA and albumins[END_REF]. With respect to above, HSA binding constant of the NiL2 is in a good range (2-10 × 10 4 ) and also is comparable to analogue compounds and some potent drugs [START_REF] Kazemi | Synthesis, characterization and biological application of four novel metal-Schiff base complexes derived from allylamine and their interactions with human serum albumin: experimental, molecular docking and ONIOM computational study[END_REF].

In addition, our compound is uncharged which is an advantage for a drug. The researches show that the unionized drugs are soluble in lipid and are capable to cross through the membrane's lipid bilayer while the ionized analogue species fail to cross [START_REF] Dighe | Synthesis of β-carboline-based n-heterocyclic carbenes and their antiproliferative and antimetastatic activities against human breast cancer cells[END_REF].

Circular dichroism (CD) spectroscopy

The CD spectra of HSA and DNA in free and bonded to nano-scale NiL2 complex states are shown in Fig. 5. The CD spectrum of HSA indicates two negative peaks at 208 and 222 nm which are characteristic absorption peaks of α-helical structure of protein attributed to the n-π* transfer for the bonds of the α-helix [START_REF] Zsila | Probing the binding of the flavonoid, quercetin to human serum albumin by circular dichroism, electronic absorption spectroscopy and molecular modeling methods[END_REF]. According to Fig. 5, the intensity of HSA peak has been altered slightly due to binding nano-complex. Thus, there is no considerable structural change for HSA due to its interaction with NiL2. The CD spectrum of free double stranded DNA exhibits four marker peaks at about 210 (negative), 222 (positive), 245 (negative) and 280 nm (positive). The positive peak at about 280 nm due to base stacking and negative peak at 245 nm due to polynucleotide helicity are main peaks of DNA for its B-type, which are highly sensitive towards the DNA interaction with small molecules [START_REF] Ivanov | Different conformations of double-stranded nucleic acid in solution as revealed by circular dichroism[END_REF]. Fig. 5 shows that the intensity of DNA peaks has been altered due to binding nano-complex more than the changes in HSA CD spectrum. This observation confirms that the interaction of NiL2 with DNA is stronger than its interaction with HSA that is in consistent with the fluorescence results.

Viscosity measurement

To further verify the interaction mode of the nano-scale NiL2 dispersion with CT-DNA, viscosity measurements of DNA upon addition of the nano-scale NiL2 dispersion were carried out. A classical intercalative mode causes an increase in the DNA viscosity due to increasing of overall length of DNA [START_REF] Kazemi | Synthesis, characterization and separation of chiral and achiral diastereomers of Schiff base Pd (II) complex: a comparative study of their DNA-and HSA-binding[END_REF]. On the other hand, groove binding and electrostatic interactions could bend the DNA helix, reduce its length and cause the reduction or no change in the DNA viscosity [START_REF] Kazemi | Synthesis, characterization and separation of chiral and achiral diastereomers of Schiff base Pd (II) complex: a comparative study of their DNA-and HSA-binding[END_REF].

The effect of the nano-scale NiL2 dispersion complex on the viscosity of CT-DNA is demonstrated in Fig. S6. As it can be seen, the viscosity of DNA decreased until molar ratio of [NiL2] / [DNA] = 0.25 and remained constant with increasing amounts of the complex, indicating the binding mode of NiL2 may be groove binding. This result is in good agreement with molecular docking and spectroscopic results.

Molecular docking studies

Molecular docking is a valuable technique to predict the stable structure of receptor-ligand complex for better recognition of the interaction details in drug discovery process. This method is frequently employed as virtual searching tools in primary steps of drug design and development.

In order to find out the preferred location of the compound on DNA and HSA, molecular docking studies were carried out. Table 3 represents the obtained binding modes and docking energies of the compound during its interaction with both of the biomacromolecules. The docked model suggests that NiL2 binds with the minor groove of DNA with standard binding free energy (ΔG°) of -7.69 kcal•mol -1 . This score is related to the docking cluster with maximum population and minimum binding energy. Also, intermolecular energy which covers H-bonding interactions, van der waals, desolvation and electrostatic energies, was estimated as -8.79 kcal•mol -1 [START_REF] Li | Studies of bicalutamide-excipients interaction by combination of molecular docking and molecular dynamics simulation[END_REF]. Also, docking studies showed that there is one H-bond between O atom of the NiL2 and DG4 nucleotide. Table 3. Molecular docking results for the interaction of NiL2 complex with DNA and HSA. The results of molecular docking for interaction of NiL2 and HSA indicate that the compound is situated in the warfarin binding site, in IIA subdomain (Table 3). The NiL2 complex is bound to the HSA with the binding energy of -5.59 kcal•mol -1 . Also, the intermolecular energy for HSA is -6.69 kcal•mol -1 . The amino acids residues with major role in the stability of NiL2-HSA system are LYS199, TRP214, ARG218, ARG222, LEU219, ALA291, ILE290, PHE223, ILE264, LEU260, ARG257, GLU153, TYR150 and ALA261. According to the obtained results, there is one hydrogen bond between N atom of the NiL2 complex with ARG222 residue of HSA. Moreover, four π-cation interactions were formed between phenyl rings of the compound with residues of the binding site. The participant residues of HSA in π-cation interactions are ARG222, LYS199 and ARG257.

Type of biomolecule binding energy (kcal•mol

Molecular dynamics simulation

The beginning structures for the MD analyses were selected from the conformations with lowest docking energies. The stability of the system (biomacromolecule and compound) properties was examined by means of RMS deviations (RMSD) of unbound and ligand conjugated biomacromolecule with respect to the initial structure, RMS fluctuations (RMSF), the solvent accessible surface area (SASA) of biomacromolecule, radius of gyration (Rg) and number of Hbonds. Preliminary simulations over 10 ns time were performed on DNA and HSA with NiL2. Table 4 provides the average values of RMSD, SASA, Rg and number of H-bonds of the HSA and HSA-NiL2 systems. Fig. 6 presents time dependence of RMSD, Rg, RMSF and H-bond values of HSA and HSA-NiL2 systems. Fig. 6A indicates that the trajectories of both systems are stable and their RMSD reached equilibrium and fluctuated around its mean value after about 5 ns simulation time. Figs. 6B and S7 show that complexation with NiL2 changed the structure of HSA when compared to free protein. It can be clearly seen that the Rg and SASA were larger upon binding of the NiL2 suggesting a less compact structure of HSA in the presence of NiL2. Decreasing the number of intramolecular H-bonds is another evidence that confirms a less compact structure of HSA in the presence of NiL2 (Fig. 6C). Also, local protein mobility was analyzed by calculating the time averaged RMSF values of free HSA and HSA-NiL2 complex and was plotted against residue numbers based on the last 5000 ps trajectory (Fig. 6D). The profiles of atomic fluctuations were found to be very similar to those of free HSA and HSA-NiL2 complex. RMSF highlights the conformational adjustments of the protein structure, in conjunction with NiL2 conformational adaptation to its binding sites. The results indicate that the residues that were in contact with the NiL2 are stable and have low RMSF values suggesting that the structure of drug binding site remains rigid during simulation. Hence, it can be concluded that the interactions of protein and NiL2 were stable during the simulation time.

QM/MM studies

In this work, a two-layer ONIOM calculations including (PM6: UFF) were employed to perform QM/MM studies. The molecular mechanics (MM) was described using the UFF force field for HSA and DNA (low layer) while semi-empirical quantum mechanics (QM) method (PM6) was chosen for the NiL2 (high layer). The most stable geometries obtained from the molecular docking calculations were opted as starting points of HSA-and DNA-NiL2 for the two-layer ONIOM studies. All calculations were performed using Gaussian 03 quantum chemistry package. The partial atomic charges on the atoms of the NiL2, DNA and HSA were considered to reoptimize the optimized geometries. The van der Waals (vdw) and electrostatic interaction were considered in structure optimization. Moreover, the polarization of the wave function of the NiL2 was considered to obtain more truthful geometries [START_REF] Kazemi | Synthesis, characterization and biological application of four novel metal-Schiff base complexes derived from allylamine and their interactions with human serum albumin: experimental, molecular docking and ONIOM computational study[END_REF]. The results of ONIOM calculations represented that structure of the NiL2 deviates from its initial geometry due to the binding to DNA and HSA. Along with the interaction of NiL2 with biomacromolecules some bond lengths and bond angles are changed. Such the changes can be indicative of the strength of the interaction between biomacromolecules and NiL2. Fig. 8, Fig. 9 show the changes take place in some geometrical parameters including bond angles of NiL2 along with binding to DNA and HSA. The calculated binding energies (ΔE) of NiL2 and biomacromolecules are present in Table 5 for DNA and HSA.

In view of considering the polarization of the NiL2 in electrostatic field of HSA and DNA in QM/MM calculation, this energy is more accurate than calculated energy by molecular docking method.

Conclusions

In conclusion, we have synthesized a new nano-scale Schiff base Ni(II) complex in water as a green solvent and at ambient temperature. The structure of complex has been established by using single crystal X-ray diffraction analysis. The crystallographic data reveal that NiL2 is placed in a crystallographic center of symmetry that is right in the middle point of the Ni(II) in a planartransoid conformation. The MMT assay results indicated that the anticancer activity of the compound is affected by its size. The DNA-and HSA-binding of the NiL2 was also investigated, using both experimental (fluorescence quenching, CD analysis, viscosity measurements and UV-Vis spectroscopy) and computational methods (molecular docking, molecular dynamics simulation and ONIOM). The obtained HSA-and DNA-binding constant values from experimental and computational methods stated that the binding of Ni(II) complex to the minor groove of DNA and warfarin binding site of protein by formation of hydrogen bond, π-cation and hydrophobic interactions. MD simulation studies revealed that complexation with NiL2 changed the structure of HSA when compared to free protein. Finally, the ONIOM results showed that the structural parameters of the compound changed along with binding to DNA and HSA, indicating the strong interaction between the compound and these biomacromolecules.

[

  Scheme 1. (A) Synthetic routes for the preparation of NiL2 in its bulk and nano-scale forms. (B) ORTEP representation of NiL2. Displacement ellipsoids are drawn at the 50% probability level and H atoms are shown as small spheres with arbitrary radii. Hydrogen bonds are shown as dashed lines.

Fig. 1 .

 1 Fig. 1. (A) Hydrodynamic size distribution of nano-scale NiL2 dispersion in water was measured by dynamic light scattering. (B) SEM image of the nano-scale NiL2 complex.

Fig. 2 .

 2 Fig. 2. The viability percentage of HeLa cancer cells after 48 h exposure to different concentrations (10, 20, 50 and 100 μM) of NiL2 in both of scales.

Fig. 3 .

 3 Fig. 3. (A) Changes in the fluorescence spectra of DNA-EB in the presence of increasing amounts of NiL2 in Tris buffer pH = 7.5 ([DNA] = 25 μM, [NiL2] = 250 μM). (B) Plot of Ln((F0 -F)F) versus Ln[NiL2]. (C) Absorption spectra of nano-scale NiL2 complex in the presence of various amounts of DNA ([NiL2] = 200 μM, [DNA] = 1 mM). (D) Plot of 1 / ɛa -ɛf versus (1 / [DNA])).

Fig. 4 .

 4 Fig. 4. (A) Changes in the fluorescence spectra of HSA in the presence of increasing amounts of nano-scale NiL2 in phosphate buffer pH = 7 ([HSA] = 10 μM, [NiL2] = 200 μM). (B) Plot of Ln((F0 -F) / F) versus Ln[NiL2]. (C) Absorption spectra of HSA in the presence of various amounts of nano-scale NiL2 ([HSA] = 10 μM, [NiL2] = 200 μM). (D) Plot of 1 / A -A0 versus 1 / [NiL2].

Fig. 5 .

 5 Fig. 5. The CD spectra of HSA and DNA in the absence and presence of nano-scale NiL2 complex in the molar ratio [complex] / [biomolecule] = 1.

- 1 )

 1 View of binding site Zoom in binding

Fig. 6 .

 6 Fig. 6. (A) Time dependence of RMSD. RMSD values for HSA and HSA-NiL2 complex during 10 ns MD simulation. (B) Time-evolution of radius of gyration (Rg) for HSA and HSA-NiL2 complex during 10 ns MD simulation. (C) The number of Intramolecular H-bonds for HSA and HSA-NiL2 complex during 10 ns MD simulation. (D) RMSF values of HSA and HSA-NiL2 complex plotted against residue numbers.

Fig. 7

 7 presents time dependence of RMSD and RMSF values of DNA and DNA-NiL2 systems. Fig.7Aindicates that the trajectories of both systems are stable and their RMSD reached equilibrium and fluctuated around its mean value after about 1 ns simulation time. Also, local DNA mobility was analyzed by calculating the time averaged RMSF values of free DNA and DNA-NiL2 complex and was plotted against atom numbers based on the last 8000 ps trajectory (Fig.7B). The results indicate that the atoms of DNA bound NiL2 are stable and have lower RMSF values than the free DNA. Hence, it can be concluded that the interactions of DNA and NiL2 were stable during the simulation time.

Fig. 7 .

 7 Fig. 7. (A) Time dependence of RMSD. RMSD values for DNA and DNA-NiL2 complex during 10 ns MD simulation. (B) RMSF values of DNA and DNA-NiL2 complex plotted against atom numbers.

Fig. 8 .

 8 Fig. 8. (A) The structure of the NiL2 complex and DNA optimized considering charged embedding of MM part in QM part in the gas phase. (B) The changes of structural parameters of NiL2 during binding to DNA.

Fig. 9 .

 9 Fig. 9. (A) The structure of NiL2 complex and HSA optimized considering charged embedding of MM part in QM part in the gas phase. (B) The changes of structural parameters of NiL2 during binding to HSA.

Temperature (°C) Kb ΔG° (kJ/mol) ΔH° (kJ/mol) ΔS° (J/mol•K)

		25	2.45 × 10 4 -25.05	76.11	339.32
	CT-DNA	35	8.37 × 10 4 -29.03	76.11	341.25
		45	1.67 × 10 5 -31.82	76.11	339.28
		25	1.55 × 10 4 -23.92	61.83	287.60
	HSA	35	5.57 × 10 4 -27.99	61.83	291.51
		45	7.36 × 10 4 -29.64	61.83	287.52

Table 4 .

 4 The average values of RMS deviations (RMSD), the solvent accessible surface area (SASA), the radius of gyration (Rg) and number of H-bonds of the HSA and HSA-NiL2 systems.

	System	RMSD (nm)	SAS (nm 2 )	Rg (nm)	Number of H-bond
	HSA	0.25 ± 0.05	298.22 ± 3.70	2.69 ± 0.02	481.50 ± 10.97
	HSA-NiL2	0.28 ± 0.04	303.69 ± 3.79	2.71 ± 0.01	462.80 ± 11.20

Table 5 .

 5 The calculated binding energies (ΔE) of NiL2 with HSA and DNA.

	System	NiL2(PM6) (kcal•mol -1 )	Biomolecule(UFF) (kcal•mol -1 )	ΔE (kcal•mol -1 )
	NiL2-HSA -2.4473	-4913.25	4957.091
	NiL2-DNA -2.4473	2823.01	-2744.14

Acknowledgment

The financial supports of the research council of University of Isfahan are gratefully acknowledged.

Appendix A. Supplementary data

CCDC 1547057 contains the supplementary crystallographic data for NiL2. Copies of the data can be obtained, free of charge, on application to CCDC, 12 Union Road, Cambridge CB21EZ, UK (Fax: +44 1223 336033 or e-mail: deposit@ccdc.cam.ac.uk).


Franc ¸ois Gauthier-Clerc

Hoel Le Capitaine

Fabien Claveau

Philippe Chevrel

Operating data of a specific Aquatic Center as a Benchmark for dynamic model learning: search for a valid prediction model over an 8-hour horizon

This article presents an identification benchmark based on data from a public swimming pool in operation. Such a system is both a complex process and easily understandable by all with regard to the stakes. Ultimately, the objective is to reduce the energy bill while maintaining the level of quality of service. This objective is general in scope and is not limited to public swimming pools. This can be done effectively through what is known as economic predictive control. This type of advanced control is based on a process model. It is the aim of this article and the considered benchmark to show that such a dynamic model can be obtained from operating data. For this, operational data is formatted and shared, and model quality indicators are proposed. On this basis, the first identification results illustrate the results obtained by a linear multivariable model on the one hand, and by a neural dynamic model on the other hand. The benchmark calls for other proposals and results from control and data scientists for comparison.

I. INTRODUCTION

The need for energy and commodities is constantly increasing and must be carefully managed to avoid possible future shortages due to limited resources. From an energetic perspective, we observe a switch from highly controllable electric plants powered by fossil fuel to renewable non-controllable sources like wind and solar. This paradigm shift will introduce more discrepancies between demand and production, leading to more volatile prices. From a physical/hardware aspect, improving existing plants, or building new modern efficient plants, become more and more expensive in the context of commodities inflation. Despite the fact that they lead to a lot of savings, these methods let existing systems aside.

These observations show the need to improve actual plant efficiency by retrofitting process control to allow more consumption flexibility and global efficiency. This approach does not involve a lot of investment since we are dealing with software upgrades only.

Today, a lot of energy-intensive plants are controlled by a simple controller without any economic optimality concern [START_REF] Drgoňa | All you need to know about model predictive control for buildings[END_REF]. Although numerous optimal controllers are available in the {firstname}.{surname}@univ-nantes.fr academic world [START_REF] Marín | Use of a predictive control to improve the energy efficiency in indoor swimming pools using solar thermal energy[END_REF], [START_REF] Alajmi | Energy audit of an educational building in a hot summer climate[END_REF], their deployments are still expensive and not well-established in some parts of the industrial world. This situation is even more pronounced in small and medium facilities like swimming pools or small HVAC (Heating Ventilation, and Air Conditioning) systems due to a lack of resources allocated to control. The challenge is therefore to propose advanced control methods that are efficient, have low development/deployment costs and are suitable for a wide variety of industrial systems, taking advantage of their common characteristics. The temperature management of swimming pools is part of this scope.

Today, many regulation systems (especially with swimming pools) are based on a closed-loop PID system with a fixed working point. It is a reliable and cheap regulation, which explains its popularity. For the swimming pool case, we observe a constant 28-degree target most of the time [START_REF]Fina facilities rules[END_REF]. Although this strategy is popular, its lack of economical optimality tends to question its use in the context of energy restrictions (implying an increasing cost) and non-permanent temperature demand (closed at night). More advanced expert based-technique are available to achieve better regulation. They rely on scheduling methods to reduce the average temperature in order to make savings (decreasing the average temperature can save up to 7% energy demand per degree [START_REF] Marín | Dynamic simulation model and empirical validation for estimating thermal energy demand in indoor swimming pools[END_REF]). They are based on static [START_REF] Escrivá-Escrivá | Application of an energy management and control system to assess the potential of different control strategies in hvac systems[END_REF] or dynamics temperature [START_REF] Lee | A data-driven method for determining zone temperature trajectories that minimize peak electrical demand[END_REF], [START_REF] Zemtsov | Economic mpc based on lpv model for thermostatically controlled loads[END_REF] scheduling methods to update the setpoint according to demand.

Unfortunately, their approaches are still simple and mainly rely on expert knowledge to formulate strategies. The validity of open loop schedulings is affected by exogenous effects like outdoor temperature or other disturbances [START_REF] Henze | Predictive optimal control of active and passive building thermal storage inventory[END_REF] which make them expensive to deploy and maintain. This situation highlights the need for other techniques, like model-based control. Instead of using expert knowledge for the control design, we can use a data-based system identification method [START_REF] Ala'raj | Data-driven based hvac optimisation approaches: A systematic literature review[END_REF] that allows model-based control methods. These methods generate a model-based control directly on the constraints (humanly interpretable). This approach avoids the need for experts and therefore reduces the deployment cost for equivalent accuracy in some tasks.

In the context of system control with economical awareness, the EMPC (Economical Model Predictive Control) algorithm has gained more and more interest in recent years. As explained earlier, this control method does not focus on setpoint (like tracking MPC) but uses an economical function instead. There are several contributions to HVAC systems [START_REF] Drgoňa | All you need to know about model predictive control for buildings[END_REF] and swimming pool [START_REF] Marín | Use of a predictive control to improve the energy efficiency in indoor swimming pools using solar thermal energy[END_REF]. For medium and small facilities like 978-1-5386-5541-2/18/$31.00 ©2018 IEEE swimming pool, this control method is simple to do (small state space) and allow economical saving [START_REF] Caspari | The integration of scheduling and control: Top-down vs. bottom-up[END_REF].

In the field of swimming pools, there is a lack of study that involves data-driven approaches to control. However since some interactions, like evaporation, are difficult to model [START_REF] Shah | Improved method for calculating evaporation from indoor water pools[END_REF], data-based methods are even more justified. Some work seems to show this interest in the context of public swimming pools [START_REF] Lovell | Thermal performance prediction of outdoor swimming pools[END_REF].

In order to help the academic community to develop new frugal energy management solutions, this article defines the outline of a benchmark, by making available the operating data of a public swimming pool. It makes it possible to consider the problem of learning dynamic models for control, MPC in particular. The proposed dataset is interesting for several reasons. It comes from the sensors which are conventionally fitted to this type of installation, in operation, over a long period. This dataset also contains different seasons which have an important impact on the energy balance of the installation and add a realistic complexity. To initiate the benchmark, we propose two identification results as a baseline for future contributions.

This document is organized as follows; Section 2 introduces the swimming pool system with a description of the use case and some overviews of the physical models, for a better understanding of the process. Section 3 details the data set and assessment criteria based on the need for data-driven control. Section 4 presents the primary results obtained for two selected black-box models.

II. THE PUBLIC SWIMMING POOL SYSTEM: WATER TEMPERATURE REGULATION

A. Benchmark case study

A swimming pool is a plant that needs to regulate pool temperature and quality as well as air humidity and temperature. It uses several local regulations to maintain water quality, air state and water pool temperature (see Fig. 1). For the air and pool temperatures, the plant uses a heating system that consists of a primary hot water loop to distribute energy across the plant from a given source of energy (mostly a dedicated boiler). Boilers are regulated thanks to some local controllers. In practice, each regulation is done with a local hardware regulator provided by each sub-part manufacturer. Global regulations for this type of facility are very rare or even non-existent.

The water constantly circulates from each pool to its recycling system in order to: heat, clean and store the water. The water quality is maintained by UV lamps and filters. The system also includes a buffer tank to store the water when people are using the pool (less water volume is needed). The water is also renewed with new clean water. The fresh water flow rate is based on the measured filter pressure and is activated 2 or 3 times a day for a short period of time. This water refilling is the main perturbation of this system, as the fresh water is very cold compared to the pool temperature. The water is heated by a countercurrent heat exchanger and the temperature is regulated according to temperature sensors before the heat exchanger. The temperature regulation is based on a 3-way valve that bypasses or not the water flow rate in the hot side of the exchanger. The water flow rate in the recycling loop is maintained with a pump and overflow mechanism. The flow rate setpoint is often constant but can change according to static scheduling (day/night). This flow rate can affect the dynamic of the heating process.

The aimed facility is a public swimming pool located in Brittany, FRANCE. This swimming pool is composed of two main pools (sport and recreational) as well as one hot tub and paddling pool. The sports pool has a volume of 672 m 3 (12.5m × 25m) with 5 lanes. This pool is coupled with a buffer tank of 30m 3 . The second pool (of 212 m 3 ) is coupled with the paddling pool (of 12m

TABLE I: Swimming pool characteristics

The swimming pool is open during the whole year from 7 a.m to 5:30 p.m each day. We observed between 250 and 600 swimmers across the day. The temperature requirement is constant during the opening period, 28 degrees for the sports pool and 31 for the recreational one. On certain Wednesdays throughout the year, baby swimming lessons are organized and require a high temperature of 32 degrees in the second pool. All specifications are available in the TABLE I.

In this facility, the heating system consists of a wood boiler and a gas boiler. The first one is used as a priority, and the second one is used when there is a shortage of wood or a high energy demand. The wood boiler provides a max power of 650kW and the gas boiler a power of 350kW. A peak energy production (with both boilers) can reach around 1MW.

In addition to the double boiler system, a sports hall uses the energy network for its air heating system. This particular energy consumption is monitored with dedicated sensors.

This public swimming pool gathers all the particularities of such a system (coupled pool, dual boiler, external energy consumption). It is complex enough to require special attention if traditional expert approaches are used to design the supervision system, regardless of the methodology considered; simulation-based approach or expert control rules. One can expect data-driven approaches to be able to handle all the details in an automatic process.

B. Identification needs for local optimal regulation

Based on our earlier motivation, the easiest and cheapest way to improve these simple regulations is to improve only one or two regulators. This approach reduces the deployment cost and can save a large part of the possible gain. From this idea, The diagram represents the water system only. The water pumps are not shown due to lack of space. In the real system, there are several radiator networks as well as small auxiliary pools using the same system as pool 1.

it is obvious to take control of the pool temperature regulation in order to save energy. We choose to take care of the two larger pools of the facility; Pool 1 (sports pool) and Pool 2. The Pool 3, the paddling pool regulation is ignored since it represents only 5% of the second pool volume.

The pool temperature can be manipulated with the control variable (3-way valve) in order to reduce the average temperature (reduce energy losses) or store energy when it is cheap while staying compliant with the user's requirements (thanks to predictive control).

The main action to save energy is to reduce the temperature when the pool is closed. In order to do this optimally, it is important to reheat the temperature early enough to achieve a compliant temperature at the pool opening. In this context, it is crucial to have a model that can predict the right time to heat or choose the pool with the highest energy-saving potential (if we have a limited boiler power). In this context, an 8-hour prediction according to the potential heating time is chosen for this benchmark. This horizon is large enough to handle the morning heating which can be a slow process. The model also needs to take into account exogenous data (such as outside temperature, air moisture, etc.) to be reliable over different periods. Moreover, the open loop 3-way valve and temperature are chosen to allow more efficient control (be able to saturate the control for instance). For an illustration, please refer to [START_REF] Caspari | The integration of scheduling and control: Top-down vs. bottom-up[END_REF].

To be able to manage all pools, the model needs to take into account the set of pools as a whole. Indeed, boilers have limited power which means that pools can interact with each other through energy consumption. This situation leads us to a MIMO formulation in order to deal with their problems and remain in a purely black box approach.

Before presenting the dataset, let us describe the theoretical model of pool energy balance which helps to understand macroscopic physical phenomena and data selection. This theoretical model, although enlightening, is insufficient. It will not be used in the process of identifying the system, the hypothesis retained being to proceed from the data; without having recourse to physical knowledge of the process to be controlled, and thus keep the cost of deployment cheap.

C. Theoretical aspect

The pool dynamics can be summarized by its heat balance, taking into account all sources of heat loss and gain.

V pool ρ w c w dT pool dt = Qtot (1)
In a traditional indoor swimming pool, we can sum up all of these interactions as follows:

Qtot = Qevap + Qrad + Qcond + Qconv + Qrefill + QHE (2)
From the equation above, we have evaporation ( Qevap ), radiation ( Qrad ), conduction ( Qcond ), convection ( Qconv ), water recycling (or refill) ( Qrefill ) and heat exchanger ( QHE ) contribution to the temperature dynamic ( see [START_REF] Calise | Energy and economic analysis of energy savings measures in a swimming pool centre by means of dynamic simulations[END_REF] for more details). Based on these interactions, the evaporation phenomenon is highly dependent on the activity of the pool, which is very difficult to measure [START_REF] Asdrubali | A scale model to evaluate water evaporation from indoor swimming pools[END_REF]. The water recycling and the heat exchanger contribution are driven by the local controller. Like water evaporation, the impact of fresh water depends on its temperature, which is rarely measured in practice and is often seasonally correlated.

The heat exchanger allows for heat transfer from the primary loop to the pool. This 3-way valve is controlled with an opening setpoint that induces a non-linear static relationship with the water flow rate [START_REF] Zajic | Modelling of an air handling unit: A hammerstein-bilinear model identification approach[END_REF]. The hot water flow rate gives energy to the pool according to the recycling flow rate. This interaction can be approximated using the NTU method [START_REF] Incropera | Fundamentals of heat and mass transfer[END_REF]. Energy transfer can be measured using upstream and downstream sensors, although depending on the system, the upstream sensor may not always be available (since conventional control does not require it). The primary water circuit that carries the heat is also connected to other subsystems such as HVAC, water tanks (for showers), and radiators (see Fig. 1). This situation can lead to a heat shortage at some point if all of those heat consumers are calling for energy. The circuit is made in parallel in order to avoid prioritizing certain subsystems.

The air state has its own energy and humidity balance. This dynamic includes the pool hall regulation system, pool state, and heat loss to the outside. This temperature is regulated with a local controller to maintain a constant setpoint.

We summarize all of the elements in Fig. 2, where the input signals are on the left side and the output signals are on the right.

III. BENCHMARK DATASET

A. Data selection

The main objective is to predict pool temperature using only available sensors to minimise deployment costs. We aim to use as few exogenous signals as possible to ensure that the identification can be deployed to a large number of installations. In this context, we include the total boiler power, [START_REF] Alajmi | Energy audit of an educational building in a hot summer climate[END_REF] Fig. 3: MIMO formulation is measured using a dedicated sensor that can encapsulate the complexity of the heating system. If this sensor is not available, it can be deduced from the input/output temperature and the flow rate of the primary loop. The air temperature and humidity are relevant for heat losses with air. The 3-way valve is the control signal, and the new water flow rate is a perturbation. Despite all this information, some data are not easy to capture and are not included in this study. For example, wind speed [START_REF] Marín | Use of a predictive control to improve the energy efficiency in indoor swimming pools using solar thermal energy[END_REF] is necessary for evaporation modeling, as well as pool attendance, but they are not available here. We also include the sports hall energy consumption in this study since it can affect the heating power for the pools. This variable is the only one that is specific to this installation. The sensor selection is summarized in Fig. 3.

The temperature of the primary loop is not included as it is part of the system state and must be used as an output signal. Including it would increase the complexity of the identification task and the integration into predictive control algorithms. This choice respects the purely black box data-driven approach and forces the identification task to include the energy transport mechanism in an implicit manner. On the other hand, the air state is more stable and can be predicted with an ad-hoc algorithm (like time series prediction [START_REF] Mendoza-Serrano | Smart grid coordination in building hvac systems: Empc and the impact of forecasting[END_REF]).

The data is subsampled using a moving average, and the signal is processed to remove sensor faults. Due to missing data and special events like COVID-19 lockdown (see section III-B), the dataset is split. The data is also normalized according to common practice [START_REF] Ljung | System Identification: Theory for the User[END_REF]. The timeline in Fig. 4 summarises all this information as well as the type of control.

The dataset includes both seasons, winter and summer. The temperature setpoint for pool 1 fluctuates between 27 and 28 degrees, while for pool 2, it fluctuates between 30 and 31 degrees. The data was collected using several hardware automatons, with each signal being measured every minute and stored since 2018. The benchmark uses data from September 2019 to August 2020. We keep this range since the swimming pool has been extended with an outside pool that complicates the analysis (more exogenous data to take into account). Moreover, this period contains several special events that generated good data diversity for model validation. Among these events, the COVID-19 period forced the pool to close for 3 months, which created a period with a low temperature. In addition, the pool uses two control strategies: one PI regulator with a constant setpoint and a more advanced strategy that reduces the temperature at night. This second method allows for better data for identification since the control signal is not correlated with the feedback signal, and it closely approximates the optimal strategy. Both types of regulation can be compared in Fig. 5. In this figure, Pool 2 is regulated according to a static setpoint, while Pool 1 is regulated with a dynamic setpoint that allows for a reduced temperature at night. The static setpoint regulation does not provide much information, as the only perturbation to the temperature is from the addition of new fresh water. In contrast, the dynamic regulation strategy adds a lot of diversity and implements the energy-saving strategy mentioned earlier.

Furthermore, this period includes four situations where the temperature was managed differently than usual, which we refer to as "scenarios". These scenarios were caused by factors such as human error and sensor faults and can be considered as potential planning scenarios for model-based implicit controllers such as EMPC. Two of these scenarios (scenario 1 and scenario 3) for the first pool are represented in Fig. 6. Note that these scenarios were excluded from the training dataset.

The four scenarios are as follows:

• A stuck 3-way-valve on pool 1 at 50%.

• Pool 2 cooled down to 25 degrees while pool 1 remained at its nominal temperature. • Pool 1 cooled down to 25 degrees while pool 2 remained at 25 degrees. • Both pool 1 and pool 2 were heated up to their nominal temperature at the same time period.

C. Identification problem & Evaluation

The identification problem is closely linked to the use that will be made by the control (MPC) of the resulting model. The identification is therefore evaluated through a prediction horizon which is not restricted to one step, as is conventionally

L(I, J) = 1 J -I + 1 J i=I K-1 k=0 ||ŷ[i + k|k] -y[i + k]|| 2 2 K
(3) The prediction ŷ[k + i|k] is estimated based on the observation at instant k (y[k]) and exogenous and control signals from the instant k to k + i. In this benchmark, the vector y contains the temperatures of the two pools. This horizon accuracy is obtained over the entire test set in a convolutive manner.

It is useful to visualize the accuracy of each candidate model according to the prediction depth by plotting the function x → L(x, x) for x in [1, H] (Fig. 7).

From the function L, we need to define several criteria to validate and evaluate the model for control. To do so, we have three sections of testing data placed at different system regimes during the year of data. Secondly, we have four scenarios to evaluate the model in an extrapolation context.

In order to use L and enrich the analysis, three criteria have been defined to evaluate the model's prediction ability for simulation and control. The first criterion includes the entire prediction window, the second is limited to short-term II).

In addition to evaluating the model's performance on the test set, analyzing the model's extrapolation capacity is also critical in this black-box identification for control. To analyze this capacity, we use the full horizon accuracy per scenario to determine which extrapolation scenario is well-predicted or not compared to the test set.

To summarize, the goal is to achieve good accuracy in the 8hour horizon with both short-term and long-term accuracy. The model's ability to extrapolate is critical to assess its usefulness in the control scheme. Prior physical knowledge can guide the identification process, but it is important to keep in mind the initial goal of cheap deployment across different swimming pool architectures.

All resources for this benchmark are available at https:// benchmark-datadriven-sysid.purecontrol.com. These resources include data and tools for plotting and evaluating models with Python and Matlab.

IV. IDENTIFICATION METHOD & PROPOSED MODELS AND RESULTS

A. State of art

System identification is a vast field that includes model estimation, model reduction, experiment design and more [START_REF] Ljung | System Identification: Theory for the User[END_REF]. In this paper, we focus on formulating a suitable model in the context of black box modeling. Unlike other methods, we mainly rely on real data, which is more challenging to deal with. Among these black box models, we can cite linear polynomial model (ARX, OE, ARMAX), linear state space (LSS) [START_REF] Ljung | System Identification: Theory for the User[END_REF], neural network (NN), linear parameter varying (LPV) [START_REF] Thilker | Identification of non-linear autoregressive models with exogenous inputs for room air temperature modelling[END_REF] and many more. This large dictionary allows users to choose the appropriate model complexity depending on the problem.

Identification of linear models has been carried out in both simulation [START_REF] Ma | Economic model predictive control for building energy systems[END_REF] and real-world contexts [START_REF] Virk | Model-based control for hvac applications[END_REF] for HVAC systems using the ARX formulation. These data-driven models are used in EMPC to manage indoor temperature based on endogenous and exogenous signals, with models capable of predicting up to 24 hours ahead. In a purely data-driven approach, the authors of [START_REF] Ferkl | Ceiling radiant cooling: Comparison of armax and subspace identification modelling methods[END_REF] utilized subspace methods to identify an entire building using sensors.

With greater expressivity, we can mention data-driven LPV (Linear Parameter Varying) identification [START_REF] Thilker | Identification of non-linear autoregressive models with exogenous inputs for room air temperature modelling[END_REF] using parametric methods (cubic spline) to model non-linear exogenous effects in an HVAC context.

In recent years, pure non-linear models like neural networks have gained more attention. They are considered universal approximators that allow for a wide range of applications.

Many applications in the control of HVAC systems can be cited [START_REF] Afram | Artificial neural network (ann) based model predictive control (mpc) and optimization of hvac systems: A state of the art review and case study of a residential hvac system[END_REF]. In this paper, the authors highlight the fact that despite the computational burden due to non-linear MPC, their application is possible with modern computers for control. Neural Networks under NLARX (Non-Linear Auto Regressive eXogenous) modeling have been successfully identified with real data [START_REF] Afram | Artificial neural network (ann) based model predictive control (mpc) and optimization of hvac systems: A state of the art review and case study of a residential hvac system[END_REF], [START_REF] Mechaqrane | A comparison of linear and neural network arx models applied to a prediction of the indoor temperature of a building[END_REF]. Neural network models often outperform linear models in terms of pure accuracy metrics [START_REF] Afram | Artificial neural network (ann) based model predictive control (mpc) and optimization of hvac systems: A state of the art review and case study of a residential hvac system[END_REF], [START_REF] Mechaqrane | A comparison of linear and neural network arx models applied to a prediction of the indoor temperature of a building[END_REF]. However, it is important to note that the generalizability property of these models with different data distributions is not guaranteed.

As a first identification application on this proposed benchmark, it seemed to us interesting to select both a pure linear state-space model with the subspace method and a neural network-based NLARX model.

There is also a large number of hybrid models in the grey box paradigm capable of achieving good results through the addition of knowledge [START_REF] Zajic | Modelling of an air handling unit: A hammerstein-bilinear model identification approach[END_REF], [START_REF] Thilker | Identification of non-linear autoregressive models with exogenous inputs for room air temperature modelling[END_REF]. However, for now, we are setting them aside as we try to avoid knowledge-based formulations as much as possible in order to stay within the paradigm of a cost-effective and generic approach in order to reduce deployment costs.

B. Identification procedure

Let us now successively consider the identification of the process by a linear dynamic model and a nonlinear NLARX model respectively.

a) Linear state space estimation : First, the subspace identification method [START_REF] Ljung | System Identification: Theory for the User[END_REF] is used to get a representative state space linear model. This identification is carried out through the Matlab toolbox [START_REF]The MathWorks, System Identification Toolbox[END_REF] (ssest function). This method uses all of the previously mentioned input signals (cf. Fig. 3). The model is as follows:

x

[k + 1] = Ax[k] + Bu[k] y[k] = Cx[k] (4)
x and u are the state and input signals respectively, and A, B, and C are matrices characterizing the state transition and the observation function. To perform prediction, we use the function forecast, which performs an ordinary least square method to estimate the initial state with respect to the size of the given past horizon (we chose 20 time steps).

b) NLARX estimation: The nonlinear model considered is defined by:

y[k + 1] = f θ (y[k], • • • y[k -n a ], u[k], • • • , u[k -n b ]) (5)
This model is learned from the data by considering the following loss function for identification, which is well suited for further design of an MPC control:

1 N N -1 n=0 1 P P -1 p=0 (ŷ[t 0 (n) + p|t 0 (n), θ] -y[t 0 (n) + p]) 2 (6)
With ŷ(t 0 (n) + p|t 0 (n), θ), the model prediction given the state at t 0 (n) and the exogenous data from t 0 (n) + p to t 0 (n).

The function f θ is the neural network involved in the NLARX formulation. The prediction is done recursively with the full Fig. 7: RMSE accuracy according to prediction horizon per pool. future prediction approach [START_REF] Mendoza-Serrano | Smart grid coordination in building hvac systems: Empc and the impact of forecasting[END_REF]. We use the solver ADAM [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] (an unconstrained non-linear solver) with the Tensorflow library using the BPTT (Back Propagation Through Time) [START_REF] Werbos | Backpropagation through time: What it does and how to do it[END_REF].

We trained both linear and non-linear models with the same MIMO formulation and procedure. A hyper-optimization procedure was used using a cross-validation method to reduce bias with this relatively small dataset. The cross-validation is done with data sections instead of regression entries to avoid overlap between samples. The testing dataset is placed at the dataset extremum to avoid any data leak. All the optimized parameters are detailed in

C. First numerical results

After training, the models are evaluated on both the test and extrapolation sets. In Fig. 8, we show the temperature predictions for both pools. The dashed red lines depict the time when the predictions for both models are readjusted with the real data. Please note that tools for generating these comparative curves are also available on the benchmark website.

We can now introduce the raw prediction accuracy of the two selected models on the test set (Fig. 7). We plotted the accuracy with respect to the prediction horizon and the pool ID. As expected, the accuracy decreases as the prediction horizon increases. The neural network model achieves higher accuracy than the linear one, especially in long-term prediction. This gap is more pronounced for the second pool, which is more challenging to predict due to higher temperatures and more intense activity of people. From the table, we can see that the NLARX model has a better accuracy compared to the linear model. However, this difference is less noticeable in short-term accuracy, which highlights the need to consider horizon accuracy for better analysis in predictive control applications. Despite the good results obtained on the test set, a significant decrease in accuracy is observed in all extrapolation scenarios for the non-linear model compared to the linear model. The linear model provides a more consistent level of accuracy between the interpolation and extrapolation regimes.

This first result shows the limits of the pure non-linear model in the extrapolation regime of the learning data. This can be problematic for the control, which will be particularly lacking in resilience in this situation. The linear model seems better adapted to ensure a good regulation of the process considered in this benchmark. However, its accuracy in the training data domain is limited and may induce suboptimal control. This result is expected but underlines the need for a model that combines the qualities of the two models.

V. CONCLUSION

The energy transition requires installations that consume less energy and are more flexible in the face of fluctuations in the energy market. In this context, operating data from a public swimming pool have been made available, to be freely and widely shared. This dataset and evaluation tools have been proposed to serve as a benchmark for learning techniques. It should make it possible to design a dynamic black box model, capable of predicting the evolution of the temperature of the pools of the public swimming pool considered. This benchmark will make it possible to compare and, as a corollary, to improve identification tools based on data. It targets patterns for control; with the idea of reducing the cost of deploying advanced control methods.

The characteristics of the target system and the associated datasets are described. The problem was formalized through an input-output representation and evaluation criteria. Matlab-Python functions are offered on the benchmark website. Associated with the proposed datasets, two dynamic models were identified. Their results are quite good in the test set, but they have some weaknesses in the extrapolation regime (especially for the non-linear model).

In future work, we plan to address the practical limitations of current nonlinear models. More precisely, we aim to develop a method in the spirit of the one presented in [START_REF] Djeumou | Neural networks with physics-informed architectures and constraints for dynamical systems modeling[END_REF]. This will involve improving the performance of the nonlinear model by providing it with linear extrapolation capabilities. We encourage the research community to propose the most promising solutions and evaluate them on the benchmark's operating data for comparison purposes.

Fig. 1 :

 1 Fig. 1: Diagram of the heating network

  Fig. 4: Dataset timeline

Fig. 5 :Fig. 6 :

 56 Fig. 5: Illustration of pool regulation of the dataset

Fig. 8 :

 8 Fig.8: Five 48 steps prediction on test data. Finally, we present the results with the quantitative criteria introduced earlier for both models in TABLE V.

  3

  ). They share the same recycling flow and a buffer tank of 35 m 3 .

	Characteristic	Value	Unit
	Pool 1 volume	672	m 3
	Pool 1 depth	1.3 to 3.0	m
	Pool 1 recycling flow rate	189	m 3 /h
	Pool 1 buffer tank size	30	m 3
	Pool 2 volume	212	m 3
	Pool 2 depth	0.45 to 1.3	m
	Pool 2 recycling flow rate	142	m 3 /h
	Pool 2 buffer tank size	35	m 3

  way valve state, indoor air temperature and humidity, outdoor air temperature, recycling speed, refill water flow rate, and outside temperature as input signals. The total boiler output

	refill valve
	water flow	downstream
	rate	temperature
	3 way
	valve
	Boiler
	Power
	Fig. 2: the pool sub system
	Pools recycling flow rate
	water refill flow rate	Pool 1 feedback
	3 way valve state	temperature
	Air temperature
	Air humidity Boiler power Sport hall heat consumption	Pool 2 feedback temperature
	Outside temperature

TABLE II :

 II Model accuracy criteria prediction, and the last one is for long-term prediction. The goal is to achieve good accuracy for an 8-hour horizon with good short-term and long-term accuracy (refer to TABLE

  TABLE III and IV.

	Parameters	Values
	N# layers	{1,2,3}
	N# unit per layer	{8,16,32,64}
	learning rate	[1e-5, 1e-2]
	batch size	[16,256]
	l2 penalization	{0.0,1e-4,1e-3}
	loss function horizon	[3, 5, 8, 15]

TABLE III

 III

	: optimized NLARX model's parameters
	Parameters	Values
	n x	[2, 3, 4]
	focus	{simulation, prediction}
	subspace horizon	{48, 24, -1}
	noise model	{None, Estimate}
	TABLE IV: optimized SS model's parameters

TABLE V .

 V

		LSS	NLARX
	full horizon acc	0.32	0.23
	Short term acc	0.11	0.10
	Long term acc	0.52	0.31
	Scenario 1 acc	0.29	0.22
	Scenario 2 acc	0.49	0.60
	Scenario 3 acc	0.88	1.24
	Scenario 4 acc	0.36	1.45

TABLE V :

 V Accuracy score for the proposed models