domain: chim sdv.tox info

text:
Valérie Forest
email: vforest@emse.fr

Jean-François Hochepied

Lara Leclerc

Adeline Trouvé

Khalil Abdelkebir

Gwendoline Sarry

Vincent Augusto

Jérémie Pourchez

Towards an alternative to nano-QSAR for nanoparticle toxicity ranking in case of small datasets

Keywords: Nanoparticles, Nanotoxicology, Ranking, Small datasets, Nanodescriptors

Statistical analysis approaches have been developed to predict the biological response to nanoparticles, especially quantitative structure-activity relationship (QSAR) models. But one major limitation remains the quantitative lack of data to build accurate models. The aim of this study was to investigate if simple alternative mathematical models could rank nanoparticles in a very binary way (i.e. toxic or not) in case of small dataset. We synthesized and characterized 25 nanoparticles from 6 metal (hydr)oxide families with particle size and shape tuning. We assessed their toxicity on RAW 264.7 cells and investigated relationships with both physicochemical and dimensional descriptors. A simple partial least square (PLS) regression analysis allowed ranking nanoparticles with respect to their toxicity, without false-negative results. But this model was not predictive due to the specific response of each family to dimensional parameters variations. A classification tree extracted the same main bulk descriptor as PLS, but interestingly showed the relevance of dimensional descriptors for the second and third node.

We thus recommend the development of family-specific models and propose the combination of these two simple methods as pre-screening tools, a compromise to bridge the gap between case-by-case studies (expensive and time-consuming) and sophisticated nano-QSAR models (not suitable for small datasets).

Introduction

During the last decades, nanotechnologies and their applications have amazingly expanded leading to a growing public concern about the toxicity of manufactured nanoparticles. The huge number and variety of engineered nanoparticles makes it impossible to test them on a case-bycase basis [START_REF] Chatterjee | Calculating the costs of nanohazard testing[END_REF]. In addition, this approach would require intensive animal testing to obtain for each nanoparticle the necessary data [START_REF] Dekkers | Towards a nanospecific approach for risk assessment[END_REF]). It appears therefore necessary to develop alternative approaches to help prioritize materials that need the most rigorous testing, at the same time prioritizing financial and human resources available in toxicological investigations [START_REF] Chatterjee | Calculating the costs of nanohazard testing[END_REF][START_REF] Puzyn | Using nano-QSAR to predict the cytotoxicity of metal oxide nanoparticles[END_REF][START_REF] Gajewicz | Advancing risk assessment of engineered nanomaterials: application of computational approaches[END_REF][START_REF] Winkler | Applying quantitative structure-activity relationship approaches to nanotoxicology: current status and future potential[END_REF].

To meet such expectations, predictive models are currently being developed. In particular, the Quantitative Structure-Activity Relationship (QSAR) model represents a promising approach.

It was initially developed in the 1960's for chemical compounds [START_REF] Hansch | Correlation of Biological Activity of Phenoxyacetic Acids with Hammett Substituent Constants and Partition Coefficients[END_REF]. As an example, QSAR approach is commonly used for risk assessment of pesticides. Indeed, due to a growing need for toxicity assessment and the increasing variety and number of products, regulatory agencies and legislations, such as REACH (Registration, Evaluation and Authorisation of Chemical, the European legislation on chemicals) encourage the use of alternatives to animal testing. Computational approaches can accelerate advances in (eco)toxicological understanding as they support the experimental data with additional in silico studies and results [START_REF] Mas | Application of chemometric methods to environmental analysis of organic pollutants: A review[END_REF][START_REF] Hamadache | Environmental Toxicity of Pesticides, and Its Modeling by QSAR Approaches[END_REF][START_REF] Villaverde | Computational Methodologies for the Risk Assessment of Pesticides in the European Union[END_REF][START_REF] Villaverde | Considerations of nano-QSAR/QSPR models for nanopesticide risk assessment within the European legislative framework[END_REF]. It was recently proposed to transpose QSAR models to nanomaterials to rapidly and cheaply screen and predict their toxicity [START_REF] Pan | Nano-QSAR modeling for predicting the cytotoxicity of metal oxide nanoparticles using novel descriptors[END_REF]). This adaptation is still at the development stage and is not without challenges [START_REF] Gajewicz | Advancing risk assessment of engineered nanomaterials: application of computational approaches[END_REF]). To continue with the example of pesticides, the environmental fate and behavior of these products formulated as nano are complex and vary from traditional formulations, arguing for the need of innovative and adapted risk assessment methodologies [START_REF] Villaverde | Considerations of nano-QSAR/QSPR models for nanopesticide risk assessment within the European legislative framework[END_REF].

But the usefulness of computational approach in the nano field is constantly rising as illustrated by some interesting results reported in the literature. In particular, Puzyn et al. [START_REF] Puzyn | Using nano-QSAR to predict the cytotoxicity of metal oxide nanoparticles[END_REF][START_REF] Kar | Periodic table-based descriptors to encode cytotoxicity profile of metal oxide nanoparticles: a mechanistic QSTR approach[END_REF] have recently applied "nano-QSAR" starting from a dataset grouping 17 metallic oxide particles all differing by their chemical composition, and concluded that their toxicity to Escherichia coli bacteria could be predicted on the basis of the chemical nature only thanks to calculated or experimental descriptors such as metal electronegativity and cation charge. Surprisingly, dimensional or shape descriptors did not play any role. Gajewicz et al. (2015) confirmed these observations in a model of mammalian cells, the HaCaT cell line, regarding the toxicity of 18 metal oxide nanoparticles. More recently, [START_REF] Pan | Nano-QSAR modeling for predicting the cytotoxicity of metal oxide nanoparticles using novel descriptors[END_REF] demonstrated that in addition to the cationic charge as shown by Puzyn et al., the individual size and aggregation size of metal oxides nanoparticles were important factors for the cytotoxicity to both E. coli and HaCaT cells.

Despite these promising results, some hurdles remain to overcome. They may be roughly classified into 3 categories: i) the choice of the most suitable nanodescriptors, ii) the choice of the most relevant biological assays; and last but not least, iii) the quality and quantity of the data. Firstly, regarding the choice of the physico-chemical descriptors, many have been proposed as suitable for nanomaterials, they can be either derived from experimental data (for instance dimensional descriptors) or calculated with various theoretical approaches [START_REF] Winkler | Applying quantitative structure-activity relationship approaches to nanotoxicology: current status and future potential[END_REF]Gajewicz et al. 2015b;[START_REF] Pan | Nano-QSAR modeling for predicting the cytotoxicity of metal oxide nanoparticles using novel descriptors[END_REF]. Often, physico-chemical descriptors are calculated, not assessed experimentally, and they usually do not follow the entire list of parameters described in the ISO recommendations (ISO/TR 13014:2012). In addition, almost always, only one sample per chemical composition is used. Secondly, a successful nano-QSAR also requires well-selected biological assays, at the same time facing the specific issues and concern about nanoparticles and providing valuable data (biological descriptors) for the QSAR model [START_REF] Fourches | Quantitative nanostructure-activity relationship modeling[END_REF][START_REF] Gajewicz | Advancing risk assessment of engineered nanomaterials: application of computational approaches[END_REF][START_REF] Dekkers | Towards a nanospecific approach for risk assessment[END_REF]. Thirdly, and most of all, the main issue of nano-QSAR lies in the quality and quantity of data. A sufficient number of samples (particles with well-defined and distinct physico-chemical descriptors) should be available to build a statistically acceptable dataset [START_REF] Fourches | Quantitative nanostructure-activity relationship modeling[END_REF][START_REF] Dekkers | Towards a nanospecific approach for risk assessment[END_REF].

Indeed, traditional QSAR approach is based on very large datasets (hundreds or even thousands substances). However, available experimental data for nanoparticles are still not enough (units or tens, in the best case) and empirical data and physico-chemical and biological characterization are not systematic when they are not inconsistent [START_REF] Gajewicz | Advancing risk assessment of engineered nanomaterials: application of computational approaches[END_REF].

Thus, although nano-QSAR models have been the topic of intensive research during the last two decades their use is actually limited, mainly because of the small size of the dataset available for modeling. We may therefore wonder if alternative mathematical models, less complex and less sophisticated cannot be sufficient to rank, although more modestly but still efficiently, nanomaterials with respect to their toxicity in case of small datasets. And indeed, recently, other approaches have gained interest such as classification trees or read-across basedapproaches (Gajewicz 2017a, b;[START_REF] Gajewicz | Decision tree models to classify nanomaterials according to the DF4nanoGrouping scheme[END_REF]. In this context, the aim of the present study was to determine if tools as simple as a simple partial least square regression or decision trees could accurately rank nanoparticles (i.e. toxic vs. non toxic) based on a small dataset. To compensate the quantitative lack of data, their quality must be optimal. We therefore synthesized 25 nanoparticles from 6 metal (hydr)oxide families we then thouroughly characterized according to the ISO standards (ISO/TR 13014:2012). We next assessed their toxicity in vitro using a model of macrophages and biological endpoints used in empirical nanotoxicology studies. This is the main originality of this work along with the synthesis, systematical and empirical physico-chemical characterization of the nanoparticles. Based on these experimental data two predictive models and a decision tree were constructed and allowed us to discuss about the relevance of these approaches.

Experimental

Nanoparticle synthesis

25 nanoparticles from 6 metal (hydr)oxide families (SiO2, TiO2, CeO2, AlOOH, ZnO, Ni(OH)2) were synthesized. In each chemical family particle size and shape were tuned to ensure an independent analysis of chemical and dimensional descriptors.

The synthesis of boehmite particles was inspired from [START_REF] Yang | Synthesis of γ-Al2O3 nanowires through a boehmite precursor route[END_REF]. Boehmite particles were synthesized by room temperature precipitation of aluminium chloride with soda-ammonia solutions followed by hydrothermal ripening at various pH values (4.5, 6 and 8.5) in autoclaves at 180°C during 24h or microwave heating (Monowave 300, Anton Paar) at 180°C during 10 minutes. The objective was to control shape from platelets to rods. All samples were washed by successive centrifugations in water, then freeze-dried. Table 1 summarises the experimental conditions. AlOOH-1 and AlOOH-3 were supposed to give the same particles but samples exhibited a different specific surface area and agglomeration level in bundles, which may be due to accidental temperature oscillation in the autoclave and/or pH shift before heating. The synthesis of zinc oxide particles (slightly elongated particles or short rods) was inspired from [START_REF] Trenque | New insights into crystallite size and cell parameters correlation for ZnO nanoparticles obtained from polyol-mediated synthesis[END_REF]. Dihydrated zinc acetate was dissolved in diethyleneglycol (0.1 M for ZnO-1, 1 M for ZnO-2) then the solutions were heated at 200°C by microwave irradiation (Monowave 300) during 3 minutes.

The synthesis of cerium dioxide particles was inspired from [START_REF] Florea | Three-Dimensional Tomographic Analyses of CeO2 Nanoparticles[END_REF]. Our protocols of precipitation of cerium(III) nitrate with soda-ammonia solution followed by hydrothermal ripening were previously described [START_REF] Forest | Impact of cerium oxide nanoparticles shape on their in vitro cellular toxicity[END_REF], they allowed shape tuning from isotropic particles to nanorods.

Synthesis of titanium dioxide nanoparticles: for TiO2-1, a 0.5M titanium butoxide and 1M triethanolamine solution was heated at 180°C during 10 min by microwave irradiation. The sample was washed by successive centrifugations in water then freeze-dried. TiO2-2 is an amorphous product resulting from titanium butoxide hydrolysis at room temperature. TiO2-3

(spheres) and TiO2-4 (rods) were commercial products from Sigma-Aldrich.

The synthesis of silica nanoparticles was inspired from Music et al. (2011). A solution of sodium silicate 26.5% was diluted 4 times in water, heated at 80°C and acidified to pH=2 by dropping a 1.8M sulfuric acid solution. For SiO2-1 heating was maintained during 45 min after acidification without pH control, for SiO2-2 heating was maintained during 90 min and pH value was maintained at 2 by sulfuric acid addition during heating. The samples were washed by successive centrifugations in water.

Nanoparticle physico-chemical characterization

Following ISO guidelines (ISO/TR 13014:2012), the following physico-chemical descriptors were considered:

-nanoparticle size: average dimensions were measured by Transmission Electron Microscopy (TEM) using a FEI TECNAI 20FST operating at 200 kV. Each particle was defined by 1 (sphere) or 2 (rods/plates) dimensions (sufficient considering that a revolution axis practically always exists): dmax and dmin.

-chemical composition: known from synthesis conditions and X-ray diffraction. Associated descriptors were x(H2O) from the raw formula M(n+)O(n/2).xH2O, the oxidation degree of the metal (cation charge) n_oxy_M, the radius of the metallic cation r_cat, and the Pauling electronegativity of the metallic element EN_M. We decided to introduce x(H2O) to make a chemical difference (in addition to solubility) between oxides and hydroxides, the other 3 descriptors were chosen from their occurrence in literature.

-solubility: we considered the solubility of the bulk materials at pH=7, from tables.

-shape: shape factor was classically defined as SF=dmax/dmin from dimensions by TEM, and we also introduced a corrected shape factor to discriminate between rods and plates, which is of course critical since rods have 2 dimensions in the nano-range when plates only have 1:

considering the length of the axis of revolution (daxis) and the dimension perpendicular to this axis (dperp), the corrected shape factor was defined as CSF=log(daxis/dperp). This way, we have different domains of CSF for each shape: for rods CSF>0, spheres CSF=0, plates CSF<0.

-the zeta potential was measured at pH=7 by dispersing nanoparticles in water by ultrasonic treatment and using a nanoZS (Malvern Instrument, UK).

-the specific surface area SSA was measured on dried powders by nitrogen gas adsorption at 77K (BET method) with an ASAP2000 instrument from Micromeritics.

-the crystalline structure was determined by X-ray diffraction using a Bruker D8 diffractometer in  geometry and Co K radiation (=1.789 Å) to check there was no minor phase in our syntheses, but was not introduced as a descriptor.

-the agglomeration/aggregation state: agglomeration is difficult to quantify and there is no associated established metrology, especially when agglomeration (weak forces between particles) and aggregation (strong forces or bridges) may both occur and evolve during biological tests. To keep it simple, we just defined 3 classes from TEM observations and stability/dispersibility when handling suspensions: Agg=0: well-dispersed particles, 1= small (submicron) or fluffy agglomerates, 2= dense or big (micrometric) agglomerates and likely aggregates.

Cellular model

The RAW 264.7 cell line derived from mice peritoneal macrophages transformed by the Abelson murine leukemia virus and was provided by ATCC Cell Biology Collection (Promochem, LGC, Molsheim, France). Cells were cultured in Dulbecco Modified Eagle Medium (DMEM) complemented with 10% of fetal calf serum and 1% of penicillinstreptomycin (called DMEMc) at 37°C under a 5% carbon dioxide humidified atmosphere.

Cells were seeded in 96-well-plates (100 000 cells in 50 µL of medium per well) and were allowed to adhere overnight. Nanoparticles were diluted in cell culture medium to reach the following final concentrations: 15, 30, 60 and 120 µg/mL. Nanoparticles were added to cells and further incubated for 24h.

Empirical toxicity testing: the lactate dehydrogenase (LDH) release

To evaluate cell membrane integrity, the cellular release in the supernatant of cytoplasmic lactate dehydrogenase was assessed using the CytoTox-96™ Homogeneous Membrane Integrity Assay (Promega, Charbonnières-les-Bains, France) according to the manufacturer's instructions. The optical density of the samples was determined using a microplate reader (Multiskan RC; Thermolabsystems, Helsinki, Finland) set to 450 nm. Three independent experiments were performed, each in quadruplicate and the activity of the released LDH was reported to that of negative control cells (incubated without nanoparticles). A positive control consisted in the maximal cellular LDH released after cells lysis.

Expression of the toxicity results: biological scores

Regarding the biological assays, three independent experiments were performed in quadruplicate. Means were calculated between the 4 doses to get an average signal and this signal was normalized to the negative control (signal obtained for the cells incubated alone, i.e.

without nanoparticles). Therefore, the basal signal obtained for the control was set to 1. As explained further, a threshold of 1.5 was selected to rank particles, i.e. a nanoparticle with a score > 1.5 was considered as cytotoxic.

Partial least square regression analysis

The physicochemical descriptors included in this study were selected either because they were recommended by the ISO/TR 13014:2012 guideline or because they have been reported in the literature to be predictive of nanoparticle toxicity. A partial least square regression (PLS) analysis was performed using R software to build and compare two models. The fisrt one included the 4 composition-related "non nano" (non dimensional) descriptors and the second one considered all (12) descriptors thus mixing chemical and dimensional descriptors. The root mean square error of prediction (RMSEP) versus the number of components was used to select the best number of components (dimensions) and the coefficients of the linear equation were validated by a leave-one-out (LOO) approach. Considering the relatively low number of samples and our main objective (to predict a binary toxic/non toxic behavior considering both physico-chemical and dimensional descriptors), we prefered the LOO approach to methods based on training and validation sets, more relevant to bigger datasets. The R² was calculated to estimate the goodness of fit of our models (i.e. the ability to reproduce the data in the training set). For internal validation we used the leave-one-out-cross-validation (LOO-CV, especially recommended for small datasets [START_REF] Gramatica | Principles of QSAR models validation: internal and external[END_REF]) and calculated the Q²LOO to evaluate the robustness and predictivity of our models. Finally, as our dataset is composed of (hydr)oxide nanoparticles, the applicability domain of the models is defined as such.

Decision tree

Decision tree classification is a classic classification method in which records are labeled and classified into different levels of a response variable according to the frequency with which a set of features appear when a certain level of the response variable takes place. The algorithm implemented in this method is known as CART (Classification And Regression Tree). The main goal of the CART algorithm is to produce a classification treea set of classifier branchesthat determines the set of rules or variable interactions capable of predicting a given outcome (response variable). The objective in the construction of these classification rules is to make possible the prediction of the occurrence of such outcomes in future scenarios using the information from the set of variables and rules described by the decision tree [START_REF] Kar | Periodic table-based descriptors to encode cytotoxicity profile of metal oxide nanoparticles: a mechanistic QSTR approach[END_REF].

Considering the small size of the dataset, we applied the leave-one-out approach for cross- Finally, in order to provide a decision aid tool to predict the toxicity of any new compound, we proposed an innovative approach to build an average tree, representing the best consensus of the n decision trees.

Results

Our approach considered, in the family of metallic oxides and hydroxides, 25 nanoparticles of 6 chemical compositions with at least 2 significantly different particles (in size and/or shape)

in each chemical family to ensure an independent analysis of chemical and dimensional descriptors. Nanoparticle samples were characterized to quantify the relevant experimental descriptors. A toxicity test, commonly used in nanotoxicology, was performed: the assessment of the LDH release from RAW 264.7 cells. The results were analyzed by a home-made ranking method. Our aim was to predict a binary behavior (probably toxic versus probably non toxic) considering both physico-chemical and dimensional descriptors using a simple mathematical model. The idea was to propose a ranking tool for rapid and cheap pre-screening purpose. The biological scores were determined for global ranking purposes with usual dose range of nanoparticles used in in vitro nanotoxicology studies and did not directly take into account the impact of the nanoparticle dose.

Following ISO guidelines (ISO/TR 13014:2012), the following physico-chemical descriptors were considered: nanoparticle size, chemical composition with associated descriptors (such as the hydration rate, the oxidation degree of the metal (cation charge), the radius of the metallic cation, and the Pauling electronegativity of the metallic element), solubility, shape, the zeta potential, the specific surface area SSA, the crystalline structure, and the agglomeration/aggregation state. Table 3 details the nanoparticle descriptors as well as the biological endpoint (score from the LDH assay) forming the dataset used in this study. Since these descriptors are not normalized we also give the order of the maximum variation of the score induced by each term: 0.15 for x(H2O), 0.80 for n_oxy_metal, 1.20 for r_cation and 0.08 for EN_metal. In other words, the cationic radius and the cation charge are the key descriptors, the introduction of x(H2O) in addition to the 3 other descriptors used as a refinement to encompass all families from metallic hydroxides to oxides has a relatively low impact, even if not to be neglected, whereas EN_metal impact is very low.

R 2 =0.718 may be thought as a bit low to claim for a model, this is obviously due to the huge dispersion of experimental scores between samples with the same chemical composition, especially CeO2 and Ni(OH)2. This dispersion has nothing to do with poor data quality but is due to morphological and dimensional effects as shown in our previous study on ceria [START_REF] Forest | Impact of cerium oxide nanoparticles shape on their in vitro cellular toxicity[END_REF]). Nevertheless, from a practical point of view, one must admit that this composition based "non dimensional" approach is sufficient to predict the non cytotoxicity of particles if we set a threshold value at score=1.5 (corresponding to the threshold above which the signal is considered significantly increased compared to the negative control, i.e. cells unexposed to nanoparticles), since no false negative results are obtained (all particles with experimental score>1.5 do have a predicted score>1.5). Only 3 particles are predicted as "cytotoxic" whereas they are not: 2 samples of CeO2 family, a material exhibiting a spectacular particle shape effect on nanotoxicity as we exhibited in a previous publication [START_REF] Forest | Impact of cerium oxide nanoparticles shape on their in vitro cellular toxicity[END_REF], and the NiO sample obtained by calcining one Ni(OH)2 sample, interesting as it shows that the correction brought by x(H2O) is not totally sufficient. The point is the dispersion of experimental scores is low with "non toxic" families as SiO2, TiO2 and AlOOH whereas it is wide with cytotoxic families as CeO2 and Ni(OH)2. Clearly, our results call for introducing more TiO2 or SiO2 samples in future studies to check this assessment.

Q 2 -LOO= 0.0 is indeed very low and indicates a priori no predictability beyond our dataset, but one must also keep in mind that our approach is not to predict accurate scores but the right quadrant.

Then we tried to build a model mixing composition and dimensional descriptors. Considering now all the chosen descriptors, after a first model we observed that 3 descriptors could be discarded without damage: dmin due to its low effect, zeta potential and solubility at pH7 due to their strong correlations to other descriptors. Considering the remaining 9 descriptors, a model with 6 components (RMSEP=0.55) gives the coefficients reported in Table 4, with R 2 =0.823. Of course, due to the relatively low number of samples, there may be overfitting. Compared to the previous model, the change of the coefficients of physico-chemical descriptors is around 20%, and among the 5 dimensional descriptors, the corrected shape factor was found not significant whereas the relative weights of the 4 other ones were relatively close.

Considering that the values quantifying the agglomeration/aggregation level were arbitrarily set to 0, 1 or 2 depending on a very rough and qualitative approach, this point calls for more attention to develop appropriate metrological tools to quantify agglomeration and aggregation.

As far as Q²-LOO is concerned, one may wonder that it is even worse than in the previous "model" and negative (!) Q²-LOO = -0.3.

One may consider that the model is poor (if existing). This is obviously true but the real reason is a key result as we explain now. As seen by comparing both figures, the corrections coming from dimensional descriptors do improve the prediction for families as Ni(OH)2, ZnO, AlOOH and SiO2 but not for CeO2 and even degrade the prediction for TiO2. The predicted cytotoxicity of NiO is significantly lowered but it remains a false positive. On the one hand, this full model has no practical benefit as compared to the simple previous model if we consider a classification based on a threshold at score=1.5. But on the other hand, it clearly shows that it does not work to mix composition and dimensional descriptors for a simple reason: the impact of dimensional descriptors does exist but is clearly family dependent. This is obvious considering for instance boehmite, where in spite of huge dimension and shape differences the experimental scores were more or less the same for all samples, whereas for ceria the shape variation induced a strong score variation. The same remarks would apply for all dimensional descriptors. At the end, the model does as best as possible to fit the results but there is no hope to be predictive this way.

As a consequence, studies considering dimensional descriptors in order to predict toxicity make sense only if restricted to a given chemical family.

At the end of this PLS study, we observe that the prediction of the right quadrant is pretty good with only the cationic radius and the cation charge, and that dimensional parameters cannot help to improve the fate of the 3 false positive because we want to deal simultaneously with all families. This is a clear invitation to switch to decision tree, very suitable for our purpose and likely to extract more homogenous subgroups. Using the leave-one-out approach for cross validation, we computed 25 decision trees and reported the metrics to assess the quality of the prediction in Table 5. For each subsample we reported the prediction result on the left out particle. The proposed method has a global accuracy of 80%, a sensitivity of 80%, a specificity of 80% and an error rate of 20%. Although the prediction approach performs well, it would benefit further tests by extending the original dataset to other compounds.

The prediction model is improved using an average tree that is built using the 25 decision trees generated before. Figure 3 depicts the final average decision tree produced considering the common features included in the whole set of subsamples used to train the model. The decision tree branches include all rules and characteristics describing the different characteristics of nanoparticles. A test on the whole dataset has an accuracy, sensitivity and specificity of 100%.

The process to construct an average decision tree involved the identification of common features in all the 25 decision trees produced from the subsamples generated. Since all 25 decision trees share the same feature (metallic cation radius < 78.75) at the root of the tree, it appeared that such feature was the most discriminatory criterion for the ranking of nanoparticles. SF criterion comes second for 72% of the generated trees, followed by dmax which is considered discriminant in 68% of the generated trees. Thus we are able to provide recommendations related to nanoparticles characteristics that should be examined first to decide on their toxicity, which are metallic cation radius, SF, and dmax.

Interestingly, the cation radius was the main descriptor according to PLS, whereas SF and dmax appeared as significant dimensional parameters in the second PLS model. So in spite of all the biases there is some consistency between both approaches and there was indeed an added value to keep them both in parallel.

Discussion

Models able to predict nanoparticle toxicity based on adequate descriptors are strongly needed.

And as a matter of fact in silico models are the subject of intensive research and different computational models for nanomaterials have been reported [START_REF] Lamon | Computational models for the assessment of manufactured nanomaterials: Development of model reporting standards and mapping of the model landscape[END_REF]). Among them, QSAR models [START_REF] Fourches | Quantitative nanostructure-activity relationship modeling[END_REF][START_REF] Puzyn | Using nano-QSAR to predict the cytotoxicity of metal oxide nanoparticles[END_REF][START_REF] Winkler | Applying quantitative structure-activity relationship approaches to nanotoxicology: current status and future potential[END_REF]Gajewicz et al. 2015b;[START_REF] Pan | Nano-QSAR modeling for predicting the cytotoxicity of metal oxide nanoparticles using novel descriptors[END_REF], read-across (Gajewicz et al. 2015a[START_REF] Gajewicz | Addressing a bottle neck for regulation of nanomaterials: quantitative read-across (Nano-QRA) algorithm for cases when only limited data is available[END_REF]Gajewicz 2017a, b), neural network [START_REF] Fjodorova | The way to cover prediction for cytotoxicity for all existing nano-sized metal oxides by using neural network method[END_REF] or decision tree [START_REF] Gajewicz | Decision tree models to classify nanomaterials according to the DF4nanoGrouping scheme[END_REF] classifications. In the present study, we could not reasonably build a QSAR model because our dataset was too small.

On the contrary, we chose PLS and decisions trees because of their simplicity and because they are easily accessible to a wider audience compared to other more complex computational methods. In this context, we showed that with a simple partial least square regression analysis it was possible to predict the cytotoxicity of a small dataset of (hydr)oxide metal nanoparticles and rank them without false-negative results considering only physico-chemical descriptors.

Furthermore, a classification tree approach confirmed the high impact of chemical descriptors.

Our results are perfectly consistent with data reported in the literature [START_REF] Puzyn | Using nano-QSAR to predict the cytotoxicity of metal oxide nanoparticles[END_REF][START_REF] Kar | Periodic table-based descriptors to encode cytotoxicity profile of metal oxide nanoparticles: a mechanistic QSTR approach[END_REF]Gajewicz et al. 2015b), that claimed a much higher impact of chemical descriptors as compared to dimensional descriptors to predict toxicity. But our approach was different as we did not exclude a priori size descriptors, we included them in our model. We evidenced there is undoubtedly a sharp impact of the particle shape with cerium dioxide [START_REF] Forest | Impact of cerium oxide nanoparticles shape on their in vitro cellular toxicity[END_REF]. Since each chemical family has its own behavior, by mixing all families the correlation with dimensional descriptors is blurred. If we focus on only one family, chemical descriptors will be constant and the impact of dimensional descriptors will be correctly evaluated. But of course that means gathering many fully characterized samples differing by their size and shape for each chemical composition of practical interest.

One may regret that we used only one cell line but it has the advantage of being an eukaryote model, poorly described in the literature related to predictive models. Indeed, most existing studies have been carried out on prokaryotic systems (i.e. Escherichia coli bacteria). Regarding eukaryotic systems information is more scarce [START_REF] Ghorbanzadeh | Modeling the Cellular Uptake of Magnetofluorescent Nanoparticles in Pancreatic Cancer Cells: A Quantitative Structure Activity Relationship Study[END_REF][START_REF] Pan | Nano-QSAR modeling for predicting the cytotoxicity of metal oxide nanoparticles using novel descriptors[END_REF].

Another originality of the present work is the use of a biological assay commonly employed in standard nanotoxicology, but not so often when trying to build predictive models where more basic indicators (i.e. EC50) are usually used, being more or less informative. The LDH release reflects cytotoxicity and it appeared to be a relevant endpoint allowing classifying nanoparticles either in a toxic or in a safe category. We have tried to construct other models including two other parameters commonly used in standard nanotoxicology: the TNF- production indicative of a pro-inflammatory response and the reactive oxygen species production signaling an oxidative stress (see Supplementary information). Suitable models could not be constructed, suggesting that the choice of the biological endpoints is not trivial and should be carefully considered, not all are relevant. Similarly, decision trees built with these data were not conclusive: when TNF- production or reactive oxygen species production were tested as target variables, key performance indicators were not satisfying (data not shown). Last but not least, another asset of our study is the use of well characterized nanoparticles, with biological assays carried out in the same conditions contrary to studies where values are extracted from previously published studies and where data obtained from different sources are directly compared. This is a major issue as it is acknowledged that a fundamental limitation for any predictive model is the quality of the data used to build it [START_REF] Gajewicz | Advancing risk assessment of engineered nanomaterials: application of computational approaches[END_REF].

And indeed, the choice of descriptors is both crucial and delicate. Experimentally measured descriptors are rarely present and all the ISO standard parameters are not systematically integrated in predictive models described in the literature while they are mandatory to consider in any standard nanotoxicology study. An inappropriate choice of descriptors can be a source of bias in the results as we can only find correlations with initially included parameters while we can miss some important finding just because the suitable descriptors were not included a priori.

Moreover, the idea of building family-specific predictive models could be a good compromise between the unsound search for a universal model able to predict the toxicity of all types of nanoparticles and the case-by-case study (that is too fastidious, highly time and moneyconsuming). Application of CART demonstrated good results for a specific target variable that is LDH release. It allowed to conclude on good predictive variables among physico-chemical characteristics of the nanoparticles. For further experiment, results provided using decisions trees can be used as a decision aid tool to decide on the most important features to measure experimentally. Also, the leave-one-out approach prevents the bias related to isolated training and test experiments and allows an objective validation of the provided decision aid tool as presented by [START_REF] Gajewicz | Decision tree models to classify nanomaterials according to the DF4nanoGrouping scheme[END_REF]. Decision tree classification appears as complementary to partial least square regression analysis as this latter allows a very binary classification of the nanoparticles (i.e. toxic or not toxic) while decision trees allow the identification of the main parameters involved in this toxicity. This could be of major interest, especially in the context of a "safer by design" approach.

Conclusion

The extension of traditional QSAR approach to nanoparticles faces some challenges, one of the biggest being the too small available datasets. With this study, based on the use of 25 thoroughly characterized metal oxide and hydroxide nanoparticles, we evidenced by partial least square regression analysis and decision tree classification relationships between toxicity and selected physical and chemical descriptors, able to rank nanoparticles in regard to their cytotoxicity without false-negative results. As illustrated by Figure 4, we thus propose these simple methods as pre-screening tools, a compromise able to bridge the gap between case-by-case studies (too fastidious, expensive and time-consuming) and more sophisticated QSAR models (efficient but not suitable for small datasets). These mathematical tools are very simple and can be easily extended to other types of nanoparticles, allowing to get new insights in the nano-(eco)toxicology field. Besides the assessment of cytotoxicity by the LDH release assay, the pro-inflammatory response and the oxidative stress were investigated through the assessment of the TNF- and reactive oxygen species (ROS) production respectively.

Figure captions

Methods

As for the LDH assay, the RAW 264.7 cells were seeded in 96-well-plates (100 000 cells in 50 µL of medium per well) and were allowed to adhere overnight. Nanoparticles were diluted in cell culture medium to reach the following final concentrations: 15, 30, 60 and 120 µg/mL.

Nanoparticles were added to cells.

TNF- production -After a 24h incubation with nanoparticles, the production of TNF-α was assessed in the supernatant using a commercial ELISA Kit (Quantikine® Mouse TNF- Immunoassay; R&D Systems, Lille, France) according to the manufacturer's instructions. The optical density of each sample was determined using a microplate reader (Multiskan RC;

Thermolabsystems, Helsinki, Finland) set to 450 nm. A standard curve was established and results were expressed in picograms of TNF- per milliliter of supernatant. Three independent experiments were performed, each in quadruplicate and the production of TNF- was reported to that of control cells (incubated without nanoparticles).

Reactive oxygen species (ROS) production -A large array of ROS activity can be assessed with the OxiSelect™ ROS Assay Kit (Euromedex, Mundolsheim, France). The assay uses the conversion of a non-fluorescent substrate, 2.7′-dichlorodihydrofluorescein diacetate that can easily diffuse through cell membranes and be converted into a fluorogenic molecule 2′.7′dichlorodihydrofluorescein (DCF) in presence of ROS: fluorescence amount is directly related to ROS level. Fluorescence was detected using a Fluoroskan Ascent fluorometer (Ex: 480 nm, Em: 530 nm, Thermolabsystems) after a 90 min and a 24h incubation of cells with the nanoparticles. A positive control was included incubating cells with H2O2 (1 mM). Three independent experiments were performed, each in quadruplicate and the generation of ROS was reported to that of the negative control (cells incubated without nanoparticles).

Expression of the results: biological scores -As for the LDH assay, means were calculated between the 4 doses to get an average signal and this signal was normalized to the negative control (signal obtained for the cells incubated alone, i.e. without nanoparticles). Therefore, the basal signal obtained for the control was set to 1.

Results

PLS analysis

Sample TNF- ROS 90min ROS 24h

AlOOH

  validation. Each learning set is created by taking all the samples except one, the test set being the sample left out. Thus, for n samples, we have n different training sets and n different tests set. This cross-validation procedure does not waste much data as only one sample is removed from the training set. We implemented all algorithms using the Python language and the SciKit learn library for machine learning methods. To assess the performance of the prediction, we provided the following metrics: accuracy, sensitivity, specificity and error rate: Accuracy = (tp + tn)/(tp+tn+fp+fn) Sensitivity = tp/(tp+fn) Specificity = tn/(tn+fp) Error rate = fn/(fn+tp) where tp, fp, tn and fn denote the number of true positives, false positives, true negatives and false negatives respectively.

Fig. 1 Fig. 2 Fig. 3 Fig. 4 Figure 1 Figure 4

 123414 Fig. 1 Plot of experimentally determined versus predicted LDH scores with model 1

Table 1 -

 1 Experimental conditions for the synthesis of boehmite particles.

	Sample	pH	Heating device	Ripening
	AlOOH-1	4.5	Autoclave	24h
	AlOOH-2	8.5	Autoclave	24h
	AlOOH-3	4.5	Autoclave	24h
	AlOOH-4	6	Autoclave	24h
	AlOOH-5	6	Microwave	10 min
	Synthesis of nickel hydroxide particles: nickel hydroxide nanoparticles are nanoplatelets, a
	shape induced by the brucite structure. For samples Ni(OH)2-1 to 4, nickel hydroxide was
	precipitated in a vortex rapid mixer (whose geometry was described by Di Patrizio et al. (2016))

Table 2 -

 2 Soda/ammonia ratios for the synthesis of nickel hydroxide particles.

	Sample	[NaOH]/[NH3]
	Ni(OH)2-1	0.1/0.3
	Ni(OH)2-2	0.1/0
	Ni(OH)2-3	0.1/0.3*
	Ni(OH)2-4	0.01/0.3
	*For sample Ni(OH)2-3, ammonia was pre-mixed with nickel nitrate solution then the solution
	was mixed with soda in the rapid mixer.

Table 3 -

 3 Dataset used for the construction of the models.

	Sample	Hydratio n rate	Oxidation degree of the metal	Radius of the metallic cation	Pauling electronegativity	Solubility	Zeta potential pH=7	Specific surface area	Dimensions	Shape factor =dmax/dmin	Corrected shape factor=log(daxis/ dperp)	Agglomeration/ aggregation state
		x (H2O) n_oxy_M	r_cat	EN_M	s(log Mtot)	zeta	SSA	dmin	dmax	SF	CSF	Agg	LDH
	AlOOH-1	0.5	3	67.5	1.61	-6.5	37	121	10	300	30	1.48		1.1
	AlOOH-2	0.5	3	67.5	1.61	-6.5	39	107	10	30	3	-0.48		1.08
	AlOOH-3	0.5	3	67.5	1.61	-6.5	37	62	6	300	50	1.7		1.05
	AlOOH-4	0.5	3	67.5	1.61	-6.5	33	40	10	100	10	1		1.07
	AlOOH-5	0.5	3	67.5	1.61	-6.5	37	54	5	50	10	1		1.07
	Ni(OH)2-1	1	2	83	1.91	-3.5	40	190	3	10	3	-0.48		2.76
	Ni(OH)2-2	1	2	83	1.91	-3.5	44	226	2	10	5	-0.7		2.57
	Ni(OH)2-3	1	2	83	1.91	-3.5	31	192	3	30	10	-1		2.99
	Ni(OH)2-4	1	2	83	1.91	-3.5	31	86	10	60	6	-0.78		2.2
	Ni(OH)2-5	1	2	83	1.91	-3.5	35	10	50	200	4	-0.6		2.07
	NiO	0	2	83	1.91	-12	35	40	10	30	3	-0.48		0.88
	ZnO-1	0	2	88	1.65	-2.5	15	93	10	10	1	0		2.62
	ZnO-2	0	2	88	1.65	-2.5	29	19	40	150	4	0.6		2.81
	CeO2-1	0	4	101	1.12	-9.5	-26	151	4	40	10	1		2.02
	CeO2-2	0	4	101	1.12	-9.5	-30	52	5	5	1	0		2.58
	CeO2-3	0	4	101	1.12	-9.5	-35	61	8	40	5	0.7		2.08
	CeO2-4	0	4	101	1.12	-9.5	-17	60	18	18	1	0		0.98
	CeO2-5	0	4	101	1.12	-9.5	-31	37	8	90	10	1		2.08
	CeO2-6	0	4	101	1.12	-9.5	9	61	15	15	1	0		0.82
	TiO2-1	0	4	74.5	1.54	-9	9	94	20	35	2	0.3		0.88
	TiO2-2	0	4	74.5	1.54	-9	1	302	2	2	1	0		1.01
	TiO2-3	0	4	74.5	1.54	-9	-4	45	20	20	1	0		1.29
	TiO2-4	0	4	74.5	1.54	-9	0	44	10	50	5	0.7		1

Table 4 -

 4 Descriptors chosen for the tentative second model and corresponding coefficients (not normalized).

	Variable	constant x_H2O n_oxy_M r_cat EN_M	SSA	dmax	SF	CSF	Agg
	coeff	-0.662	0.190	-0.485	0.037 0.104 0.00193 0.00248 -0.0111 -0.0234	0.344
	variation range		0.19	0.97	2.11	0.08	0.84	0.74	0.54	0.06	0.79
	of the term
	variable*coeff

Figure 2 plots the experimental versus the predicted LDH score with this second model.

Table 5 -

 5 Performance of LDH prediction. TP: true positive, TN: true negative, FP: false positive, FN: false negative.

	Sub	TP	TN	FP	FN
	sample
	1		1
	2		1
	3		1
	4		1
	5		1
	6			1
	7	1
	8	1
	9	1
	10	1
	11				1
	12			1
	13	1
	14	1
	15			1
	16	1
	17				1
	18	1
	19		1
	20		1
	21		1
	22		1
	23		1
	24		1
	25		1
	Total	8	12	3	2
	Accuracy	0.8
	Sensitivity	0.8
	Specificity	0.8
	Error rate	0.2

funding from Institut Mines Telecom, ROQSTAR grant.

Conflict of Interest

The authors declare that they have no conflict of interest.


Konstantin Avrachenkov
email: konstantin.avrachenkov@inria.fr

Jasper Goseling
email: j.goseling@utwente.nl

Berksan Serbetci
email: serbetci@eurecom.fr

Distributed Cooperative Caching for Utility Maximization of VoD Systems

We consider caching of VoD contents in a cellular network in which each base station is equipped with a cache. Videos are partitioned into chunks according to a layered coding mechanism and the goal is to place chunks in caches such that the expected utility is maximized. The utility depends on the quality at which a user is requesting a file and the chunks that are available. We impose alpha-fairness across files and qualities. We develop a distributed asynchronous algorithm for deciding which chunks to store in which cache.

I. INTRODUCTION

A promising means to increase the efficiency of Video-on-Demand (VoD) services in a cellular network is to proactively cache data in the base stations: part of the popular content will be stored at the base stations and the backhaul will be used only when the stored content is refreshed preferably at offpeak hours. The challenge is then to optimally place content in the base stations.

In this paper, we look at the problem from a utility perspective and develop a low-complexity distributed and asynchronous content placement algorithm for VoD systems. More specifically, we consider a layered coding (LC) mechanism, enabling to serve videos at different qualities. Users request content at a specified quality and the resulting utility is inversely proportional to the number of bits that need to be downloaded over the backhaul due to non-cached layers for the requested quality. The rationale is that latency, and therefore, quality of experience, is proportional to the size of this download. We impose alpha-fairness across files and qualities [START_REF] Srikant | The mathematics of internet congestion control[END_REF]. We provide an algorithm that optimizes the content placement in caches while keeping the communication between caches in the network to a minimum as the caches need to exchange information only with their neighbours, i.e., with those whom they have overlapping coverage area with.

Next, we provide a brief discussion work that is most closely related to the current paper. In [START_REF] Dehghan | A utility optimization approach to network cache design[END_REF], Dehghan et al. have provided a utility-driven caching for a single cache, where each content has been associated with a utility, which is a function of the corresponding content hit probability. In [START_REF] Avrachenkov | A low-complexity approach to distributed cooperative caching with geographic constraints[END_REF] we have developed a distributed asynchronous algorithm that deals with the miss probability minimization by casting the problem into the framework of potential games for general networks, and showed that our algorithm converges to a Nash equilibrium, and in fact to the best Nash equilibrium in most practical scenarios. In [START_REF] Avrachenkov | Distributed cooperative caching for VoD with geographic constraints[END_REF] the developed algorithm has been used for the residual bandwidth minimization of VoD systems with several video partitioning mechanisms. In [START_REF] Applegate | Optimal content Placement for a large-scale VoD system[END_REF], Applegate et al. have presented an approach for optimal content placement for a large-scale VoD by formulating the problem as a mixed integer problem. The main difference between the present work; and [START_REF] Dehghan | A utility optimization approach to network cache design[END_REF] is the consideration of a network with many caches instead of a single cache, and [START_REF] Avrachenkov | A low-complexity approach to distributed cooperative caching with geographic constraints[END_REF], [START_REF] Avrachenkov | Distributed cooperative caching for VoD with geographic constraints[END_REF] is the considered performance measure, and [START_REF] Applegate | Optimal content Placement for a large-scale VoD system[END_REF] is the concept of video partitioning and the heterogeneity of the video popularities over the network. See [START_REF] Avrachenkov | A low-complexity approach to distributed cooperative caching with geographic constraints[END_REF], [START_REF] Paschos | The role of caching in future communication systems and networks[END_REF] for a more exhaustive review of the existing caching work.

Our main contributions in this paper are as follows:

• We use layered coding (LC) mechanism to divide VoD contents into chunks and consider alpha-fair utility functions to represent a versatile notion of fairness for user satisfaction throughout the whole network; • We provide a distributed asynchronous algorithm which can be interpreted as giving the best response dynamics in a potential game. • We provide the optimal solution to the best response dynamics via Lagrangian; • We evaluate our algorithm through numerical examples.

We study the optimal content placement probabilities for different degrees of fairness. Let us outline the organization of the paper. In Section II we give the formal model and problem definitions. In Section III we provide the game formulation of the problem, analyze the structures of the best response dynamics and Nash equilibria. In Section IV, we present practical implementations of our low-complexity algorithm and show the resulting optimal content placements for different notions of fairness of utility in a real network.

II. MODEL AND PROBLEM DEFINITION

We consider a network of N base stations that are located in the plane R 2 . We will use the notation [1 : N ] = {1, . . . , N } and Θ = P ([1 : N ]) \ ∅, where P ([1 : N ]) is the power set of [1 : N ]. We specify the geometric configuration of the network through A s , s ∈ Θ, which denotes the area of the plane that is covered only by the caches in subset s, namely A s = (∩ ∈s Ā ) ∩ (∩ ∈s Āc ), where Ā is the complete coverage region of cache .

As a special case we will consider the case that all base stations have the same circular coverage region with radius r. In this case we specify the location of each base station, with x m for the location of base station m ∈ [1 : N ]. We then obtain Ām as the disc of radius r around x m .

Each base station is equipped with a cache that can be used to store videos from a content library C v = {c 1 , . . . , c J }, where J < ∞. Each element c j represents a video, where j is the file (video) index. The size of video c j is denoted by w j .

Next, videos are partitioned into chunks. We assume that video c j is partitioned into Q chunks and consequently, the chunk library consists of video chunks C c = {c 1,1 , . . . , c 1,Q , . . . , c J,1 , . . . , c J,Q }, where J, Q < ∞. Each element c j,q represents a video chunk, where j is the file (video) index and q is the chunk index. The chunk size for the element c j,q is denoted by w j,q .

The motivation behind partitioning videos into chunks is as follows. In this work we will consider LC mechanism [START_REF] Wiegand | Overview of the H.264/AVC video coding standard[END_REF]. If the video is partitioned by using LC mechanisms (e.g., MPEG-2, MPEG-4, etc.), the video chunks represent either the base layer or one of the enhancement layers. The base layer is necessary for the media stream to be decoded. Accordingly, the further enhancement layers are applied to improve the video quality. Hence, for the layered coding mechanism, the element c j,1 represents the base layer for video j, and c j,2 is the first enhancement layer required to obtain a better video quality, and so forth. For the LC mechanism, the chunk size w j,q will then depend on the average base size of the different video qualities, and the total size of the video j with the highest video quality is equal to

w j = Q q=1 w j,q . ( 1
)
Our interest is in users located in the plane over the area that is covered by the base stations, i.e., uniformly distributed in A cov = ∪ s∈Θ A s . The probability of a user in the plane being covered by caches s ∈ Θ (and is not covered by additional caches) is denoted by

p s = |A s |/|A cov |. A user located in A s , s ∈ Θ can connect
to all caches in subset s, and has access to all video chunks stored in these caches.

We assume that there are R video quality levels. The probability that video j is requested with video quality ρ by a user at x u ∈ s is denoted by a j,ρ,s . Here ρ = 1 refers to the video with the lowest quality and ρ = R refers to the case where the video has the highest quality. There is an obvious relation between the required chunks and the desired video quality. For the LC mechanism, files can be divided into Q = R layers with the proper size selections (w j,q ) and receiving all layers gives the video with the highest quality.

We assume that the distribution for the requested video's quality depends on the user's location, and the requested video quality for a user located in s ∈ Θ follows a known and fixed probability mass function (pmf) f (s) (ρ).

For the file popularities, we assume that a 1 ≥ a 2 ≥ • • • ≥ a J . We assume that video popularities do not vary as ρ changes. Even though any popularity distribution can be used, most of our numerical results will be based on the Zipf distribution for the video popularities for any video quality ρ. In that case, for any video quality indicator ρ, the probability that a user will ask for content j is equal to

a j = j -γ J j=1 j -γ , (2)
where γ > 0 is the Zipf parameter. Since a j and f (s) (ρ) represent statistically independent random variables, we conclude that the probability that video j is requested by a user located at x u ∈ s with video quality ρ is equal to

a j,ρ,s = j -γ J j=1 j -γ f (s) (ρ). (3)
Content is placed in caches using knowledge of the request statistics a j,ρ,s , but without knowing the actual request made by the user. We denote the placement policy for cache m as b (m)

j,q := 1, if c j,q is stored in cache m, 0, if c j,q is not stored in cache m, (4)
and the overall placement strategy for cache m as

B (m) =     b (m) 1,1 . . . b (m) J,1 . . . . . . . . . b (m) 1,Q . . . b (m) J,Q    
as a Q × J matrix. The overall placement strategy for the network is denoted by B = B (1) ; . . . ; B (N ) as an Q×J ×N three-dimensional matrix. Caches have capacity K, i.e.,

J j=1 Q q=1 w j,q b (m) j,q ≤ K, ∀m.
For clarity of presentation, we assume homogeneous capacity for the caches. However, our work can immediately be extended to the network topologies where caches have different capacities (i.e., for the case where cache m ∈ [1 : N ] has capacity K m ). Next,we are interested in designing a cache placement strategy that optimizes the sum of utilities over all files throughout the network for LC mechanism.

We will use α-fair utility functions, each yielding different notions of fairness [START_REF] Srikant | Communication networks: An optimization, control, and stochastic networks perspective[END_REF]. Then, the sum of utilities over all files throughout the network is given by

U (B) = s∈Θ J j=1 Q ρ=1 p s a j,ρ,s h 1-αρ,s j,ρ,s 1 -α ρ,s , (5)
where α ρ,s is the fairness parameter for the videos requested with quality ρ by the users located in subregion s ∈ Θ, and

h j,ρ,s = ρ q=1 w j,q 1 - ∈s 1 -b ( ) j,q ρ q=1 w j,q . ( 6
)
Here h j,ρ,s defines the portion of the requested video j with video quality ρ that is available to the user located in s ∈ Θ. Since the data available at the cache do not give rise to any load in the backhaul, the transmission rate has a direct relation with this parameter and the content stored in caches will be transmitted to the user with lower latency.

Our goal is to find the optimal placement strategy maximizing the sum of utilities of the contents throughout the overall network as follows:

Problem 1. max U (B) s.t. J j=1 Q q=1 w j,q b (m) j,q ≤ K, ∀m, (7) b
(m) j,q ∈ [0, 1], ∀j, q, m. (8)
It is easy to verify that Problem 1 is not concave. We will provide a distributed asynchronous algorithm to address Problem 1 in which we iteratively update the placement policy at each cache. In [START_REF] Avrachenkov | A low-complexity approach to distributed cooperative caching with geographic constraints[END_REF], [START_REF] Avrachenkov | Distributed cooperative caching for VoD with geographic constraints[END_REF] we have showed that we can define an algorithm that can be viewed as the best response dynamics in a potential game for similar non-concave problems. We will present a similar algorithm here. We make use of the following notation. Denote by B (-m) the placement policies of all caches except cache m. We will write U (B (m) , B (-m) ) to denote U (B). Also, for the sake of simplicity for the potential game formulation that will be presented in the following section, let U (m) denote the sum of utilities of the users within coverage region of cache m, i.e.,

U (m) B (m) , B (-m) = s∈Θ m∈s J j=1 R ρ=1 p s a j,ρ,s h 1-αρ,s j,ρ,s 1 -α ρ,s = s∈Θ m∈s J j=1 R ρ=1 p s a j,ρ,s ρ q=1 wj,q 1-1-b (m) j,q ζ (m) (j,q) ρ q=1 wj,q 1-αρ,s 1 -α ρ,s , = J j=1 Q q=1 s∈Θ m∈s p s a j,q,s q t=1 wj,t 1-1-b (m) j,t ζ (m) (j,t) q t=1 wj,t 1-αρ,s 1 -α ρ,s = J j=1 Q q=1 s∈Θ m∈s p s a j,q,s η 1-αρ,s j,q,s 1 -α ρ,s (9)
where

η j,q,s = q t=1 w j,t 1 -1 -b (m) j,t ζ (m) (j, t) q t=1 w j,t , (10)
and

ζ (m) (j, q) = ∈s\{m} 1 -b ( ) j,q . (11)

III. POTENTIAL GAME FORMULATION

In this section we provide a distributed asynchronous algorithm to address Problem 1 in which we iteratively update the placement policy at each cache. We will show that this algorithm can be formulated as providing the best response dynamics in a potential game.

In our algorithm, each cache tries selfishly to maximize the utilities of the users within its coverage region U (m) defined in [START_REF] Monderer | Potential games[END_REF]. Given a placement B (-m) by the other caches, cache m solves for

B (m) in Problem 2. max U (m) B (m) , B (-m) s.t. J j=1 Q q=1 w j,q b (m) j,q ≤ K, (12) b
(m) j,q ∈ [0, 1], ∀j, q. ( 13
)
Each cache continues to optimize its placement strategy until no further improvements can be made. At this point B is a Nash equilibrium strategy that maximizes the overall utility satisfying

U (m) (B (m) , B (-m) ) ≥ U (m) ( B(m) , B(-m) ), ∀m, B (m) .
(14) We will refer to this game as the content placement game and demonstrate in the next subsections that this game is a potential game [START_REF] Monderer | Potential games[END_REF] with many nice properties.

A. Convergence analysis

In this section, we will prove for that if the caches repeatedly update their placement strategies it is guaranteed to converge to a Nash equilibrium in finite time. Note that the order of the updates is not important as long as all caches are scheduled infinitely often.

Theorem 1. The content placement game defined by payoff function ( 9) is a potential game with the potential function given in [START_REF] Błaszczyszyn | Optimal geographic caching in cellular networks[END_REF]. If we schedule each cache infinitely often, the best response dynamics converges to a Nash equilibrium in finite time.

Proof. In order to show that the game is a potential with the potential function U (B), it is easy to verify that

U (m) B(m) , B (-m) -U (m) B (m) , B (-m) = U B(m) , B (-m) -U B (m) , B (-m) ,
which completes the proof of the first statement. The proof states that the improvement in the utility by the best response after each update is equal to the improvement in the utility in the overall network. The detailed analysis is trivial and skipped due to space constraints. Now, since there exists only a finite number of placement strategies, none of the caches will be missed in the long-run. Moreover, each non-trivial best response provides a positive improvement in the potential function in a potential game. Hence, we are guaranteed to converge to a Nash equilibrium in finite time for both games.

B. Structure of the best response dynamics

In this subsection we will analyze the structure of the best response dynamics. We will show that solution to Problem 2 can be obtained by solving a concave optimization problem and we provide the general solution for different fairness parameters α ρ,s via Lagrangian duality.

Since the feasible solution set to Problem 2 is convex and the objective function is strictly concave and continuous, the optimal solution exists and this optimal solution is a unique maximizer of the function. We can write the Lagrangian function as

L(B, λ, M, N) = U (m) B (m) , B (-m) -λ   J j=1 Q q=1 w j,q b (m) j,q -K   + J j=1 Q q=1 µ j,q b (m) j,q - J j=1 Q q=1 ν j,q b (m) j,q -1 , where B (m) , M, N ∈ R J×Q +
, and λ ∈ R. In order to achieve the maximum in L(B, λ, M, N), B (m) must satisfy

∂U (m) B (m) , B (-m) ∂b (m) j,q
-λw j,q + µ j,q -ν j,q = 0, [START_REF] Uzun | Semantic Modeling and Enrichment of Mobile and WiFi Network Data[END_REF] which is easy to solve since the Lagrangian is separable with respect to b (m) j,q . From the nature of the problem, b (m) j,q 's can take values between 0 and 1; which yields to a probabilistic placement in caches. This yields to the conclusion that one way of possibly storing chunks in the caches is by using the obtained optimal probabilities based on the optimal solution at the off-peak hours [START_REF] Avrachenkov | Optimization of caching devices with geometric constraints[END_REF], [START_REF] Błaszczyszyn | Optimal geographic caching in cellular networks[END_REF], [START_REF] Serbetci | Optimal geographical caching in heterogeneous cellular networks[END_REF]. Another strategy is to use time-to-live (TTL) caches [START_REF] Fofack | Performance evaluation of hierarchical TTL-based cache networks[END_REF], [START_REF] Dehghan | A utility optimization approach to network cache design[END_REF] by setting timers to control the storage probabilities of different chunks.

C. Structure of Nash equilibria

In this subsection we provide insight into the structure of the Nash equilibria of the content placement game. We know from the previous subsection that the games is a potential game. Hence, the Nash equilibria for the game corresponds to the optimal placement strategies satisfy the solutions of the dual problem of Problem 1.

Corollary 1. Let B denote a placement strategy at a Nash equilibrium of the content placement game. Then, B satisfies the solution obtained via [START_REF] Uzun | Semantic Modeling and Enrichment of Mobile and WiFi Network Data[END_REF] for any α ρ,s and ∀m = 1, . . . , N .

D. A real wireless network: Berlin network

In this section we will evaluate our theoretical results for the topology of a real wireless network. We have taken the positions of the base stations provided by the OpenMobileNetwork project [START_REF] Uzun | Semantic Modeling and Enrichment of Mobile and WiFi Network Data[END_REF]. The base stations are located in the area 974×812 ms around the TU-Berlin campus. We will consider these base stations as our caches with certain cache capacities. The coverage radius of the base stations is equal to r = 700 m. The locations of the base stations and their corresponding coverage areas are shown in Figure 1.

IV. PERFORMANCE EVALUATION

In this section we will present practical implementations of our algorithm for the content placement game and evaluate our theoretical results according to a network of caches with their geographical locations following a real wireless network.

A. The ROBR algorithm

We will use the Random Order Best Response (ROBR) algorithm [START_REF] Avrachenkov | A low-complexity approach to distributed cooperative caching with geographic constraints[END_REF] for the content placement game to maximize the user utilities. The basic idea of our algorithm is to repeatedly perform best response dynamics presented in Section III-B. ,s = 0 q = 1 q = 2 q = 3 q = 4 q = 5 0 50 100 150 200

File index (j) 10 -2 10 -1 10 0

,s = 1 q = 1 q = 2 q = 3 q = 4 q = 5 0 50 100 150 200

File index (j) 10 -2 10 -1 10 0

,s = 2 q = 1 q = 2 q = 3 q = 4 q = 5 0 50 100 150 200

File index (j) ,s q = 1 q = 2 q = 3 q = 4 q = 5 2. The cached portions of the chunks for different notions of fairness. For ROBR algorithm, at each iteration step, a random cache is chosen uniformly from the set [1 : N ] and updated by applying best response dynamics. We assume that all caches are initially empty (i.e., B (m) = 0 Q,J , where 0

Q,J represents a Q × J zero matrix, ∀m ∈ [1 : N ].
). The algorithm stops when U (m) B (m) , B (-m) converges. ROBR algorithm is shown in Algorithm 1.

We consider the content library of size J = 200. We assume a Zipf distribution for the video popularities, setting γ = 1 and taking a j according to [START_REF] Avrachenkov | Optimization of caching devices with geometric constraints[END_REF]. We consider the case where videos have R = 5 different video qualities. We use the video quality bandwidth requirements and the base size data given in [START_REF] Krishnappa | DASHing YouTube: An analysis of using DASH in YouTube video service[END_REF] for the chunk sizes, where the base size is the average of the sizes of a large set of tracked videos in certain video qualities given in megabytes per minute (MB/min). We assume that videos in the content library are all 40 minutes long. The corresponding video quality specifications are shown in Table I. The first column indicates the video quality index ρ, the second column indicates the base sizes of the videos and the third column indicates the sizes of the 40 minutes long videos.

The videos are partitioned into Q = 5 chunks by using LC mechanism. From the third column of Table I, it immediately follows that w j,1 = 102.4 MB, w j,2 = 69.6 MB, w j,3 = 99.6 MB, w j,4 = 348.0 MB, and w j,5 = 688.4 MB, ∀j ∈ [1, J]. Finally, we set K = 6.54 GB, i.e., each cache can store five 1080p movies.

Throughout the network we use a uniform distribution for the requested video qualities, i.e., f (s) (ρ) = 1/5, ρ ∈ 1, 2, 3, 4, 5.

In Figure 2 the cached portions of the chunks throughout the whole network is shown. For α ρ,s = 0, the optimal placement strategy is designed in order to minimize the residual bandwidth. Therefore, both sizes of the chunks and the file popularities are affecting the optimal placement. As α ρ,s increases, the effect of video popularities on the optimal solution becomes less effective; in fact, when α ρ,s → ∞, the video popularities becomes completely ineffective and the chunks for different videos at the same layer are all stored with equal probability, yielding max-min fairness. Note that in this case the chunk sizes are still taken into account since the chunk sizes are already embedded in h j,ρ,s .

Algorithm 1 :Fig. 1 .

 11 Fig. 1. Locations of Base Stations from OpenMobileNetwork dataset.