Mehdi Sahihi
email: mehdi.sahihi@gmail.com

Fatma Gaci

Isabelle Navizet

Identification of new alpha-synuclein fibrillogenesis inhibitor using in silico structure-based virtual screening

Keywords: Parkinson disease, a-Synuclein, Dihydromyricetin, Virtual screening, MD simulation

Identification of new alpha-synuclein fibrillogenesis inhibitor using in silico structure-based virtual screening

Mehdi Sahihi, Fatma Gaci, Isabelle Navizet

Highlights:

• Similar compounds to DHM were selected from PubChem library using VS method.

• αSN monomer binding affinity of the filtered molecules were analyzed in details.

• The best compound was selected as αSN aggregation Inhibitor for further analysis.

• MD simulation study showed that the inhibitor weakens the structure of αSN trimer.

Introduction

One of the age-related neurodegenerative disorders is Parkinson's disease (PD) of which the frequency increases with the increase in age. About 2-3% of people over 65 years old worldwide suffer from the disease [START_REF] Poewe | Parkinson disease[END_REF]. This disease characterized by occurrence of Lewy bodies (LBs) which is formed by abnormal aggregation and accumulation of alpha-synuclein (αSN) in existing neurons [START_REF] Lewy | Paralysis agitans. I. Pathologische anatomie[END_REF][START_REF] Xu | Alpha-synuclein in Parkinson's disease: from pathogenetic dysfunction to potential clinical application[END_REF]. αSN is a 140 amino acids protein (molecular weight of about 14 kDa) with significant conformational flexibility that can adopt a wide range of structural conformations such as fibrils and oligomers [START_REF] Deleersnijder | The remarkable conformational plasticity of alpha-synuclein: blessing or curse?[END_REF][START_REF] Mehra | α-Synuclein misfolding and aggregation: Implications in Parkinson's disease pathogenesis[END_REF]. The native structure of αSN is not still exactly defined and has been described as intrinsically disordered protein (IDP) [START_REF] Breydo | α-Synuclein misfolding and Parkinson's disease[END_REF][START_REF] Theillet | Structural disorder of monomeric α-synuclein persists in mammalian cells[END_REF], helical [START_REF] Bartels | α-Synuclein occurs physiologically as a helically folded tetramer that resists aggregation[END_REF][START_REF] Wang | A soluble α-synuclein construct forms a dynamic tetramer[END_REF], or a mixture of both [START_REF] Burré | Properties of native brain α-synuclein[END_REF]. While the function of αSN is unclear, it is highly enriched in neurons, where it is present in the presynaptic cleft as well as mitochondrial membranes. As for Alzheimer's disease (AD), for which certain oligomers of the amyloidogenic protein are toxic, the initiation of formation of αSN aggregate is suspected to be a key step in PD disease formation [START_REF] Winner | In vivo demonstration that α-synuclein oligomers are toxic[END_REF]. The cytotoxic effect of these aggregates could be because of their effect on lipid membrane permeabilisation, mitochondrial damage and oxidative stress [START_REF] Ma | Non-polyphenolic natural inhibitors of amyloid aggregation[END_REF]. Also, recent experimental studies have shown that the injection of αSN aggregates in animal's brain can stimulate their prion-like spread in other brain regions [START_REF] Sangwan | Inhibition of synucleinopathic seeding by rationally designed inhibitors[END_REF]. Most of the available drugs for PD are chemicals acting as a prodrug that will be converted in dopamine in the brain [START_REF] Chakraborty | A review of possible therapies for Parkinson's disease[END_REF][START_REF] Haddad | Dopamine and levodopa prodrugs for the treatment of Parkinson's disease[END_REF] or are dopamine agonists inhibiting the enzymes responsible of dopamine metabolism [16,[START_REF] Carrera | Current drugs and potential future neuroprotective compounds for Parkinson's disease[END_REF]. Both of these sorts of medicines cause unwanted side effects such as headaches, nausea, somnolence and so on [START_REF] Carrera | Current drugs and potential future neuroprotective compounds for Parkinson's disease[END_REF]. Hence, inhibition of αSN fibrillogenesis and disaggregation of αSN fibrils by therapeutic medications could provide a more significant neuroprotective effect against PD diseases [START_REF] El-Agnaf | A strategy for designing inhibitors of α-synuclein aggregation and toxicity as a novel treatment for Parkinson's disease and related disorders[END_REF][START_REF] Nagai | Inhibition of polyglutamine protein aggregation and cell death by novel peptides identified by phage display screening[END_REF]. Nowadays, finding of compounds that can work as robust inhibitors of αSN fibrillogenesis remains an important challenge. Given the enormous set of existing compounds, it is impossible to perform experimental screening of all available databases in a reasonable timeframe. Advances in computational chemistry provide new opportunities to steer experiments with realistic simulations on how the available compounds may be manipulated to introduce new effective inhibitor for αSN aggregation at minimum cost.

Virtual screening (VS) is known as one of the effective methods to recognize novel bioactive compounds and is being used in many drug discovery and development strategies at the moment.

In this method, the chemical entities of big libraries are rapidly screened to achieve structures with better ligand-receptor interaction pattern. In most of the VS strategies docking methods are applied on the large databases [START_REF] Vittorio | Rational design of small molecules able to inhibit α-synuclein amyloid aggregation for the treatment of Parkinson's disease[END_REF][START_REF] Jafaripour | In silico drug repositioning of FDA-approved drugs to predict new inhibitors for alpha-synuclein aggregation[END_REF].

Flavonoids are one of the large class of natural compounds that are mainly found in plant-derived materials and display antioxidant [START_REF] Van Acker | Flavonoids can replace α-tocopherol as an antioxidant[END_REF][START_REF] Yu | Antioxidant activity of citrus limonoids, flavonoids, and coumarins[END_REF] and anti-carcinogenic activities [START_REF] Kawaii | Quantitation of flavonoid constituents in citrus fruits[END_REF][START_REF] Pszczola | Natural colors: pigments of imagination[END_REF][START_REF] Cooray | [END_REF]. The ability of flavonoids to bind to carrier proteins [START_REF] Khosravi | Computational Studies on the Interaction of Arctiin and Liquiritin With β-lactoglobulin[END_REF][START_REF] Sahihi | In-silico study on the interaction of saffron ligands and beta-lactoglobulin by molecular dynamics and molecular docking approach[END_REF][START_REF] Sahihi | Binding of biguanides to β-lactoglobulin: molecular-docking and molecular dynamics simulation studies[END_REF] as well as their effects against cancer [START_REF] Williams | Flavonoids: antioxidants or signalling molecules? Free radical biology and medicine[END_REF] and nicotine-induced oxidative stress [START_REF] Sengupta | Differential roles of 3-Hydroxyflavone and 7-Hydroxyflavone against nicotine-induced oxidative stress in rat renal proximal tubule cells[END_REF] have been validated by a large number of investigations. Dihydromyricetin (DHM; Figure S1) is a flavonoid compound isolated from Ampelopsis grossedentata (a traditional southern Chinese herb), that has several physiological and pharmacological effects [START_REF] Murakami | Hepatoprotective activity of tocha, the stems and leaves of Ampelopsis grossedentata, and ampelopsin[END_REF]. Especially, DHM showed neuroprotective against oxidative stress and valuable effects in neurodegenerative diseases [START_REF] Kou | Ampelopsin attenuates 6-OHDA-induced neurotoxicity by regulating GSK-3β/NRF2/ARE signalling[END_REF]. Recently, Jia et. al. found that DHM could inhibit αSN fibrillogenesis and destabilize mature αSN fibrils in a dose-dependent manner [START_REF] Jia | Dihydromyricetin Inhibits α-Synuclein Aggregation, Disrupts Preformed Fibrils, and Protects Neuronal Cells in Culture against Amyloid-Induced Cytotoxicity[END_REF].

Their findings suggested that DHM has a great potential to be developed into a new aggregation inhibitor for αSN. Since DHM inhibits αSN fibrillogenesis and destabilizes mature αSN fibrils, its derivatives would be suitable candidates for designing of novel antiPD drugs with inhibitory effects on aggregation of αSN.

Therefore, in the present work the following computational studies have been done to design new αSN aggregation inhibitor using VS strategy:

1-More than 8600 compounds with similar chemical structure to DHM were chosen from PubChem free database [START_REF] Kim | PubChem 2019 update: improved access to chemical data[END_REF].

2-These compounds have been filtered on basis of Lipinski's "Rule of five" and topological polar surface area (TPSA) 3-The binding properties of the 314 remained compounds with αSN monomer have been investigated using molecular docking tools. 4-By considering the results of the VS part, the compound with lowest binding energy was chosen to perform a 100 ns molecular dynamics (MD) simulation of its complex with a αSN trimer to investigate its effect as potential inhibitor for αSN aggregation.

Computational procedures

2.1.Building a library of the drug-like compounds

Structure based VS is one of the common techniques being used in primary phases of drug design, discovery and development at pharmaceutical companies and academic groups [START_REF] Mirza | Discovery of selective dengue virus inhibitors using combination of molecular fingerprint-based virtual screening protocols, structure-based pharmacophore model development, molecular dynamics simulations and in vitro studies[END_REF]. To build a library of candidate molecules, the structure of DHM was selected as template and similar structures were saved from PubChem free database. The physicochemical characteristics of small molecular drugs have an important effect on their Blood-Brain barrier (BBB) permeability. Based on Lipinski's rule of five, a small molecular drug must have less than 5 H-bond donors, less than 10 H-bond acceptors, be smaller than 500 Da of molecular weight, and its calculated octanol-water partition coefficient (log P) should be less than 5 to have an appropriate BBB permeability [START_REF] Lipinski | Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings[END_REF].

The log P value of a compound is the logarithm of its partition coefficient between n-octanol and water (log (coctanol/cwater)) and is an extent of the compound's hydrophilicity. The ideal range for BBB permeability has been found to be 1.5-2.5 [START_REF] Mikitsh | Pathways for small molecule delivery to the central nervous system across the blood-brain barrier[END_REF]. We used the Lipinski's rule of five from the filters available in the PubChem database to select the candidate molecules for the VS. Then, the topological polar surface area (TPSA; as a commonly used medicinal chemistry metric for the optimization of a drug's ability to permeate cells) was calculated using Molinspiration webserver (www.molinspiration.com). Therefore, only compounds with TPSA between 20 and 130 Å 2 were kept in the VS database [START_REF] Jin | Structure-based virtual screening of influenza virus RNA polymerase inhibitors from natural compounds: molecular dynamics simulation and MM-GBSA calculation[END_REF]. DHM and 314 other ligands met the criteria.

2.2.Molecular Docking

The experimental structure of αSN saved from PDB (PDB ID: 1XQ8) database [START_REF] Ulmer | Structure and dynamics of micelle-bound human α-synuclein[END_REF]. This structure has been obtained by average structure criteria using solution NMR method. The water molecules of the .pdb file were removed and hydrogen atoms, charges, standard bond orders and the missing residues were added using Build/Check/Repair model wizard in WHAT IF web interface [START_REF] Vriend | WHAT IF: a molecular modeling and drug design program[END_REF].

Minimization was carried out up to 500 steps by steepest descent method using GROMACS 2018 package [START_REF] Berendsen | GROMACS: a message-passing parallel molecular dynamics implementation[END_REF] with GROMOS96 53A6 force field [START_REF] Oostenbrink | A biomolecular force field based on the free enthalpy of hydration and solvation: the GROMOS force-field parameter sets 53A5 and 53A6[END_REF]. The structures of the obtained drug-like compounds were saved in .sdf format and converted to .pdbqt format using PyRx 0.8 software as input files for VS [START_REF] Dallakyan | Small-molecule library screening by docking with PyRx[END_REF]. The chemical structure of DHM and the 314 other ligands (obtained after filtration by Lipinski`s rule of five and TPSA calculations) were geometry optimized by AM1 semi empirical level of theory using Gaussian 09 quantum chemistry package [START_REF] Frisch | EM64L-Gaussian 09[END_REF] and are represented in Table S1. The blind molecular docking of DHM to αSN monomer was carried out using AutoDock Vina [START_REF] Trott | AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading[END_REF] which is an academic-free molecular docking (and virtual screening) software and works based on empirical scoring functions. The molecular docking was done based on a semiflexible docking mechanism with a rigid receptor and flexible ligand. The default parameters were assigned for molecular docking calculations and the pose with the lowest binding energy was chosen as final result. For the other compounds obtained from PubChem database, the center of the docking box was set in the center of the DHM binding site and the same procedure was done.

2.3.MD simulation

The initial structure of αSN trimer used in the MD simulation has been obtained from its available amyloid fibril structure determined by cryo-electron microscopy (PDB ID: 6A6B [START_REF] Li | Amyloid fibril structure of α-synuclein determined by cryo-electron microscopy[END_REF]). All MD simulations were performed using the GROMACS 2018 package [START_REF] Berendsen | GROMACS: a message-passing parallel molecular dynamics implementation[END_REF] with GROMOS96 53A6 force field [START_REF] Oostenbrink | A biomolecular force field based on the free enthalpy of hydration and solvation: the GROMOS force-field parameter sets 53A5 and 53A6[END_REF]. The force field parameters of the selected ligand (the small molecule with the most negative binding energy to αSN monomer) were obtained from prodrg web server [START_REF] Schüttelkopf | PRODRG: a tool for high-throughput crystallography of protein-ligand complexes[END_REF] and Lemkul's method [START_REF] Lemkul | Practical considerations for building GROMOS-compatible small-molecule topologies[END_REF] was used to correct the atomic partial charges. The αSN trimer was put in a cubic box with periodic boundary conditions (PBC) in three directions. The protein trimer was placed in the center of the box and the minimum distance between solute surface and the box was set to 1.0 nm. Ten molecules of the selected ligand were randomly placed around the αSN trimer.

Then, the box was filled with SPC water molecules [50] and the system was neutralized by replacing six water molecules with Chloride ions (Cl -). In parallel, a system of αSN trimer without the ligand molecules was also prepared. The systems were energy minimized using 50000 steps of the steepest descent method and were equilibrated for 1 ns in NVT followed by 1 ns in NPT ensembles. Finally, a 100 ns MD simulation was carried out at 1 bar and 300 K for each system.

Berendsen thermostat [START_REF] Lemak | On the Berendsen thermostat[END_REF] Parrinello-Rahman barostat [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] and 1 nm cutoff for van der Walls (vdW) and Coulomb interactions were applied. Also, the long-range electrostatics were considered by particle mesh Ewald (PME) method [START_REF] Darden | Particle mesh 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] and the equations of motion were integrated by leap-frog algorithm with a time step of 2 fs. Finally, the atomic coordinates were recorded to the trajectory file every 0.5 ps for later analysis.

Results and Discussion

3.1.Virtual screening and molecular docking

Considering the structural similarity with DHM as a representative inhibitor for αSN aggregation, 8609 compounds were obtained from PubChem database. The obtained compounds were filtered using Lipinski's rule of five considering their Hydrogen bonding capacity, molecular weight, and calculated log P as well as good TPSA. Finally, 314 out of 8609 compounds were selected as potential satisfied drug with possibility of good absorption or permeability through the BBB.

The previously published experimental results suggest that IDPs remain disordered in their bound state [START_REF] Liu | Advantages of proteins being disordered[END_REF]. Hence, the αSN monomer was chosen to investigate its affinity for DHM binding. A molecular docking study was carried out to find out the preferred location and conformation of DHM on αSN monomer. The docking pose of the most stable DHM-αSN complex is represented in Figure 1. The results reveal that DHM binds to the terminal part of αSN monomer`s IDP domain which is the αSN binding site for rosmarinic acid [57]. The standard binding free energy (ΔG 0 ) of this interaction is about -5.8 kcal.mol -1 . The molecular complex system shows six H-bond interactions between DHM and Glu126, Met127, Pro128 and Glu130 residues of αSN. These Hbond interactions accompanied by hydrophobic interactions may have key roles in the stability of the system. Then, 314 compounds that passed the Lipinski`s rule of five and TPSA filtering were subjected to molecular docking (Table S1). These molecules were docked to the DHM binding site of αSN monomer. Achieved docked poses for each compound were ranked by their binding energies and conformations with the most negative binding energy were selected as main binding pose. Results revealed that 61 compounds among them showed lower binding energies than DHM. These compounds were highlighted in Table S1. The compound with PubChem ID 100968625 (IUPAC name: 5,7-dihydroxy-2-[(2R,3R)-3-(hydroxymethyl)-2-(4-hydroxyphenyl)-2,3-dihydro-1,4benzodioxin-6-yl]chromen-4-one, Figure S2) displayed the lowest binding free energy (ΔG 0 = -7.1 kcal.mol -1 ) and was selected for further analysis. The value of calculated binding free energy is very close to the obtained values for polyphenols of Ocimum sanctum (OS) extract which its inhibitory effect on aSN aggregation has been studied, before [57].

Moreover, the calculated log P value of this molecule is about 2.4 that makes it as a good BBB permeable drug. Figure 2 represents the docking pose of the most stable complex of the selected ligand and αSN monomer. Analysis of the docking results reveals that there are three H-bond interactions between the ligand and Glu126, Pro128 and Ser129 residues of αSN. All binding energies and key residues involved in the interactions of DHM, candidate molecule and the other first 21 screened compounds with the highest binding energies, with αSN monomer are summarized in Table 1. These results revealed that H-bond, hydrophobic and π-π interactions have dominant roles in binding of candidate molecules to αSN monomer.

3.2.MD simulation

The interactions between ten molecules of the selected ligand (PubChem ID: 100968625) and a αSN trimer was investigated using MD simulations method. The selected ligand (Figure S2) contains a carbonyl group, two ether oxygen atoms, and five hydroxyl groups that can interact with αSN via H-bond interaction. This makes a total of 8 possible H-bond sites for each ligand molecule. Therefore, the total number of H-bonds formed between the αSN trimer and the ten ligand molecules was calculated and shown in Figure 3. Analysis of the results showed that Hbonds were formed between the αSN trimer and ligands, rapidly after few ns of MD simulation trajectory and their number increased to about 20 at around 30 ns and fluctuated around this mean value until the end (about 2 H-bonds per ligand molecule). Figure 3b shows the representative structure of αSN and one of the ligands molecules interacting by H-bond interactions. It shows that there was a H-bond interaction between hydrogen atom of the hydroxyl group of the ligand and an oxygen atom of the carboxyl group of the side chain of the residue Glu57. Same interaction has been reported for the interaction of DHM molecules with aSN trimer [START_REF] Jia | Dihydromyricetin Inhibits α-Synuclein Aggregation, Disrupts Preformed Fibrils, and Protects Neuronal Cells in Culture against Amyloid-Induced Cytotoxicity[END_REF]. At the same time, there was another hydrogen bond between the same ligand and the amide group of residue Glu120. This observation is completely in agreement with the docking results and proves the key role of H-bond interactions in the stability of inhibitor-aSN complex. Furthermore, the previously published results have confirmed that H-bond interactions have important role in the interaction of their suggested inhibitors and aSN protein [START_REF] Jafaripour | In silico drug repositioning of FDA-approved drugs to predict new inhibitors for alpha-synuclein aggregation[END_REF][START_REF] Jia | Dihydromyricetin Inhibits α-Synuclein Aggregation, Disrupts Preformed Fibrils, and Protects Neuronal Cells in Culture against Amyloid-Induced Cytotoxicity[END_REF]57,58]. The compactness of αSN trimer during the MD simulations (in the presence and absence of the selected compound) is an important factor to study the inhibitory effect of the candidate drug and conformational changes of the protein trimer during its interaction with it. Therefore, the radius of gyration (Rg) of the aSN trimer was Inhibitor evaluated as an index of protein structural compactness. Figure 4 shows that the Rg of αSN trimer in the presence of the selected compound increased in comparison to the case without the compound indicating the less compact structure of αSN trimer (Figure 4). The average values of Rg are about 2.03 and 2.11 nm for free aSN and its complex with 10 molecules of the candidate drug, respectively. This observation proves that the compound destabilizes the compact structure of αSN trimer. To evaluate the mobility of the residues of αSN trimer, the root mean square fluctuations (RMSF) of the residues of protein aggregate in the absence and presence of inhibitor were analyzed and presented in Figure 5. The results of the RMSF analysis indicate that all the residues of αSN trimer show higher fluctuations in the presence of inhibitor molecules in comparison to the state without inhibitors.

This observation is in consistency with the Rg results and could be another evidence about good inhibitory effect of the selected compound as new inhibitor of αSN aggregation. The middle chain (chain II in figure 5) of the αSN trimer has been surrounded between two other chains and is not accessible for the inhibitor molecules as much as two other chains. This protection causes its less fluctuations in comparison to the chains I and III of the protein trimer in both states (in the presence and absence of the inhibitors).

Figure 6 represents that during 100 ns of MD simulation trajectory, most of the inhibitor molecules associate in N-terminal of the protein trimer that causes higher fluctuation of this part in three chains of αSN trimer as it is clear in Figure 5. In the presence of the inhibitor, the N-terminal fluctuations of the three chains are about 14%, 15% and 18% higher than C-terminal for Chains I to III, respectively. However, the maximum RMSF difference between the N-and C-terminal in the absence of the inhibitor would be only about 7% (for chain I). In other words, the interaction of inhibitor molecules with aSN chains weakens the intermolecular interactions between the Nterminal residues of the protein trimer and increase the fluctuations of the residues in this part of the protein. Furthermore, one inhibitor molecule inserts into the interspace between two chains of the αSN trimer, which would disturb the intermolecular interactions between the chains. To prove this claim, the distance between the Cα atoms of Glu57 of chain I and Glu120 of chain II was calculated as shown in Figure 7. The average distance in the presence of inhibitor was about 12 Å after 20 ns, while the distance in the αSN trimer in water was nearly 5 Å during the entire simulation time. Therefore, insertion of the inhibitor molecule inside the inner space of the αSN chains caused a loss of contact between two chains and disrupted their initial ordered structure.

Conclusion

In this study, we built a library of compounds with similar structure of Dihydromyricetin (DHM; a natural flavonoid with inhibitory effect on αSN fibrillogenesis). Lipinski's "Rule of five" and TPSA screenings were used to select the drug like compounds. Moreover, binding free energy and binding modes of the filtered compounds were investigated using in silico virtual screening method. The analysis of the molecular docking calculations showed that the compound with PubChem ID of 100968625 has the highest affinity for binding to αSN monomer (ΔG 0 binding = -7.1 kcal.mol -1 ) among the other candidates. Its high binding affinity for aSN suggested that this compound could be considered as an inhibitor for aSN aggregation with higher efficiency than DHM and some of the other previously suggested antiPD drugs. The behavior of aSN trimer in the presence of the selected candidate was further analyzed using molecular dynamics (MD) simulation method. Analysis of MD trajectories showed that H-bond interactions between the hydroxyl groups of the ligand molecules with side chains of the residues in αSN trimer could destabilize the compact structure of αSN trimer. On the other hand, most of the inhibitor molecules associate in N-terminal of the protein trimer and cause 100% increasing of fluctuations for the residues in N-terminal of the trimer structure. Also, investigation of Rg of the aSN trimer showed that the inhibitor molecules decrease the compactness of protein structure for about 4% in comparison to the protein free state in water. Furthermore, penetration of one inhibitor molecule into the interspace of two protein chains, disturbed the ordered structure of the aSN filament, weakened the intermolecular interactions between the protein chains and increased their distance.

In conclusion, our study revealed the inhibitory effect of selected ligand (PubChem ID=100968625) on αSN fibrillogenesis. Eventually, we believe that this study will shed the lights on finding of potentially active compounds against the αSN aggregation. We also confirm that it is necessary to do further in vitro and in vivo investigations for developing potent drug inhibitors of αSN fibrillogenesis for the treatment of PD and similar diseases.

Figure 1 .

 1 Figure 1. Docking pose of αSN-DHM complex. H-bond and hydrophobic interactions are represented as green dashed lines and red curves, respectively.

Figure 2 .

 2 Figure 2. Docking pose of the most stable complex of the selected ligand and αSN monomer. H-bond and hydrophobic interactions are represented as green dashed lines and red curves, respectively.

Figure 3 .

 3 Figure 3. Interaction of the selected ligand (PubChem ID: 100968625) with the αSN timer. (a) Total number of hydrogen bonds between the αSN trimer and the ten ligand molecules during 100 ns of MD simulation. (b) Binding conformation between residues in the αSN trimer and one of the ligand molecules at the end of trajectory. The blue solid lines represent the H-bonds. Glu residues have been shown as pink ball & stick representation.

Figure 4 .

 4 Figure 4. Time dependence of the radius of gyration (Rg) for αSN trimer during the simulation, in the absence and presence of inhibitor.

Figure 5 .

 5 Figure 5. RMSF values of residues of αSN trimer during the 100 ns of simulation time, in the absence and presence of the inhibitor. (Chain I: residues 1-63; Chain II: residues 64-126 and Chain III: residues 127-189).

Figure 6 .

 6 Figure 6. Configuration of the ligand molecules around αSN trimer during 100 ns of MD simulation.

Figure 7 .

 7 Figure 7. Change of Glu57-Glu120 distance in chains I and II of αSN trimer during the simulation, in the absence and presence of inhibitor.

[

  58]. Jayaraj, R.L. and N. Elangovan, In silico identification of potent inhibitors of alphasynuclein aggregation and its in vivo evaluation using MPTP induced Parkinson mice model. Biomedicine & Aging Pathology, 2014. 4(2): p. 147-152.

Figure S2 .

 S2 Figure S2. Chemical structure of the selected ligand from PubChem free database using virtual screening method (PubChem ID of 100968625; IUPAC name: 5,7-dihydroxy-2-[(2R,3R)-3-(hydroxymethyl)-2-(4-hydroxyphenyl)-2,3-dihydro-1,4-benzodioxin-6-yl]chromen-4-one).

Table 1 .

 1 Binding energies and key residues involved in the interactions of DHM (and the first 22 screened compounds with the highest binding energies) and αSN monomer.

Table S1 .

 S1 PubChem chemical IDs, optimized structures and binding energies of 314 drug-like chemical compounds obtained from PubChem free database. 61 compounds with lower binding energies in comparison with DHM have been highlighted.

	PubChem

ID Optimized structure Binding Energy/kcal.mol -1 PubChem ID Optimized structure Binding Energy/kcal.mol -1 TPSA/Å 2 TPSA/Å 2

	1880	-4.8	26034	-5.3
		70.67		109.36
	88281	-5.0	159026	-5.6
		90.89		74.23
	190366	-6.2	245874	-4.6
		116.84		74.24
	378687	-5.5	442410	-5.0
		68.91		86.99
	580166	-6.3	639314	-6.0
		107.60		107.60
	688853	-6.1	764011	-6.0
		90.89		96.22

Acknowledgments

M.S. is grateful to the French Agence Nationale de la Recherche (Grant: ANR-BIOLUM ANR-16-CE29-0013).

Supporting Information

Identification of new alpha-synuclein fibrillogenesis inhibitor using in silico structure-based virtual screening Mehdi Sahihi


Justine Luy

Dr Dorine Ameline

Dr Christine Thobie-Gautier

Dr Mohammed Boujtita

Estelle Lebègue
email: estelle.lebegue@univ-nantes.fr

] J Luy

Detection of Bacterial Rhamnolipid Toxin by Redox Liposome Single Impact Electrochemistry

Keywords: bacterial toxins sensing, electrochemical single impacts, redox liposomes, lipid membrane, ultramicroelectrode

The detection of Rhamnolipid virulence factor produced by Pseudomonas aeruginosa involved in nosocomial infections is reported by using the redox liposome single impact electrochemistry. Redox liposomes based on 1,2-dimyristoyl-sn-glycero-3phosphocholine as a pure phospholipid and potassium ferrocyanide as an encapsulated redox content are designed for using the interaction of the target toxin with the lipid membrane as a sensing strategy. The electrochemical sensing principle is based on the weakening of the liposomes lipid membrane upon interaction with Rhamnolipid toxin which leads upon impact at an ultramicroelectrode to the breakdown of the liposomes and the release/electrolysis of its encapsulated redox probe. We present as a proof of concept the sensitive and fast sensing of a submicromolar concentration of Rhamnolipid which is detected after less than 30 minutes of incubation with the liposomes, by the appearing of current spikes in the chronoamperometry measurement.

Pathogenic bacteria secrete various virulence factors as their biochemical weapons to kill the target host cells. They directly interact with the outer lipid bilayer membrane of eukaryotic cells, leading to the premature cell death by apoptosis or necrosis. Such virulence factors are highly involved in the bacterial infection processes, hence identifying them early is the key for the success in bacterial detection purposes. [1,2] Pseudomonas aeruginosa secretes glycolipidic biosurfactant named Rhamnolipid (RL) which have shown a broader cytolytic activity against human macrophages, keratinocytes and lung epithelial cells. [3][4][5] Therefore assays identifying RL toxin at a very early stage are of great importance for the specific detection of Pseudomonas aeruginosa which is involved in nosocomial infections. [4,6] The amphiphilic nature of RL biosurfactant allows to partition into biological and artificial membranes altering their biophysical properties. [6] Based on the RL interaction properties with the lipid membrane, we designed redox 1,2-dimyristoyl-sn-glycero-3phosphocholine (DMPC) liposomes containing 0.5 M potassium ferrocyanide as an encapsulated electrochemical active probe. Redox DMPC liposomes allowed to detect a submicromolar concentration of RL in aqueous solution in less than 30 minutes thanks to the electrochemistry of single impacts at a Pt ultramicroelectrode (UME).

The single-entity impact method is a very sensitive and powerful electroanalytical tool for biological and environmental purposes such as electrochemical detection of pathogens. [7][8][9][10] The versatility of ultramicroelectrodes (UMEs) gives unique and essential properties for electroanalytical sensing owing to their small size, high sensitivity, fast steady-state response, low double-layer charging current and minimal ohmic loss. [11] The observation of these stochastic events can provide information on various single nanoparticles contrary to ensemble (bulk) measurements. [START_REF] Bard | Nanoelectrochemistry[END_REF][START_REF] Crooks | [END_REF][14][15] The main advantage of studying collisions of single entities is the low limit of detection (in principle, one single species) inherent to this electroanalytical method and the ability to study single entities (cells, viruses, nanoparticles, etc.…) in real time (dynamic measurement). [16][START_REF] Dick | Proc. Natl. Acad. Sci[END_REF][START_REF] Dick | [END_REF][19][20][21][22][23][24][25][26][27] Electrochemistry of single redox liposome impacts at an UME consists in detecting the electrolysis of a redox probe encapsulated inside a liposome when it is released at the UME after impact (or collision). The UME is polarized at the oxidation or reduction potential of the encapsulated redox probe and a concentration of several pico-molar of redox liposomes added to the aqueous buffer electrolyte is enough to detect "current spikes" in the chronoamperometry (i-t) curve corresponding to discrete collision events. [23,[28][29][30][31][32][33][34][35][36] The electron transfer does not readily occur across a lipid bilayer, thus the electrolysis of the liposome redox active content after collision and membrane rupture or opening at the UME surface led to many studies dealing with the membrane permeation mechanism. [23,28,[30][31][32][33][34][35][36][37][38][39][40] Here, our work is based on the previous results where no current spike was observed in the chronoamperometry curve because the redox DMPC liposomes did not break during impact (or collision) onto the UME surface. [23,28] Hence, the electrochemical sensing principle is based on the weakening of the liposomes lipid membrane upon interaction with a destructive bacterial virulence factor which leads upon impact at an UME to the breakdown of the liposomes and the release/electrolysis of its encapsulated redox probe, as previously reported. [1,6,41] In the presence of RL toxin in solution (acting like a surfactant in the lipid membrane), current spikes corresponding to the electrolysis of the encapsulated redox probe released from weakened liposomes are detected (Figure 1). Thanks to the redox liposome single impact electrochemistry, the highest detection limit of RL toxin (500 nM) has been reached in comparison to several micromoles per liter previously reported. [1] Chronoamperometry measurements were recorded at +0.7 V vs Ag/AgCl (3 M NaCl) onto a 10 µm diameter disk Pt UME in 0.1 M phosphate buffer (PB) at pH 7.4 as aqueous electrolyte in the absence and presence of commercial RL (R90). The chosen concentration of RL (500 nM) presented in this work has been evaluated being the most appropriate for the proof of concept. Indeed, the toxin's detection with the redox liposomes single impact electrochemistry has to follow strict parameters such as a short time of detection (less than 1 hour after inoculation with liposomes) in order to propose a fast and sensitive sensing method.

The applied potential for the chronoamperometry measurements presented in Figure 2 is +0.7 V vs Ag/AgCl, corresponding to the steady-state current of K4Fe(CN)6 encapsulated redox probe at a Pt UME determined by cyclic voltammetry (Figure S1). The Pt UME was gently polished using a wetted fine grid silicon carbide paper from Struers (4000-grid SiC) before each electrochemical experiment. The variation of the background current (baseline) observed in the different recorded i-t curves (Figure 2) could be due to the partial covering of the Pt UME surface by liposomes and/or lipids and/or toxins. As expected, no current spike is detected in Figure 2A in all i-t curves recorded in the absence of RL, confirming that redox DMPC liposomes do not break during the impact onto Pt UME and no redox probe is released and electrolyzed. [23,28] In Figure 2A, only an increase of the background current noise is observed, probably related to the redox DMPC liposome single impacts onto Pt UME surface, in agreement with the electrochemical collision experiments by the "blocking effect" method (Figure S2). Also, in the absence of redox DMPC liposomes, in 500 nM RL aqueous solution (Figure 2B, black curve), no current spike is observed over a 30 s polarization time. The same behaviour (no current spike) is reported for the i-t curve recorded immediately after addition of an appropriate concentration of the redox DMPC liposomes solution (Figure 2B, blue curve). This suggests that an inoculation time is necessary for the interaction of RL with the liposomes lipid bilayer before to observe current spikes related to the single impact events, which is in agreement with the previous results reported in literature. [1] Indeed, the time required for observing the Rhamnolipid's effect which weak the liposomes in order to they break and/or release their content during the impact onto UME is not well understood at this stage of the study and probably dependent on both liposome and toxin concentrations. Significant current spikes appear in the i-t curve after only 15 min of incubation in RL solution (Figure 2B, orange curve), suggesting that RL affect the redox DMPC liposomes membrane leading to the breakdown during the collision onto Pt UME surface, as previously observed with Triton X-100 surfactant. [23,28] The shape of current spikes is quite broad for single collisions of 300 nm diameter liposomes containing 0.5 M ferrocyanide whose the electrolysis is fast. This is probably due to the interaction of RL toxin with the liposome lipid bilayer which leads to a slow release of the liposome redox content during the collision at the UME surface. After 30 min of redox DMPC liposomes incubation in the 500 nM RL aqueous solution (Figure 2B green curve), the current spikes are no longer clearly visible and after 50 min (Figure 2B, pink curve), the i-t curve became flat again as initially. A possible explanation of this effect is that the inoculation time is became too much high, leading to the destruction of all redox DMPC liposomes in solution in the presence of RL, hence no possible impact. Therefore the ideal time of incubation in our collisions experiment is between 10 and 30 min after addition of redox liposomes in RL solution where the current spikes are clearly distinguished in the i-t curves (Figures 2B andS3, orange curves). As a control experiments, chronoamperometry measurements were performed in the same experimental conditions but using non-redox DMPC liposomes (with only encapsulated phosphate buffer) and no current spike was observed in all i-t curves recorded over 1 hour (Figure S4). In order to confirm that the current spikes observed in the orange i-t curve (after 15 min of incubation) presented in Figure 2B are related to the redox liposome single impacts, the charge passed during the individual electrolysis was determined by integration of each current spike (Figure S5). To take into consideration that only current spikes of at least three times the background current noise (600 ± 100 fA) were considered for integration. From the measured charge of each current spike, the liposome hydrodynamic diameter was calculated using Faraday's law, as indicated in Equation 1 (Figure S5).

𝑑 lip = 2 √ 3𝑄 4𝜋𝑛 e 𝐹𝐶 redox 3

Equation 1

In Equation 1, dlip is the liposomes hydrodynamic diameter, Q is the measured charge, ne is the number of electrons transferred during the electrolysis, ne = 1, F is Faraday's constant, and Credox is the concentration of redox species encapsulated in the liposomes, Credox = 0.5 M. Here, we assume that the initial concentration (0.5 M) of encapsulated ferrocyanide is the same in each liposome and also that all of the contents is oxidized at the UME surface during the liposome collision. Based on the calculated hydrodynamic diameters, the size distribution evaluated from the electrochemical collision experiments was compared to the one obtained from dynamic light scattering (DLS) data (Figure 3). According to the overlay presented in Figure 3, the size distribution calculated from electrochemical data (orange bar) is in good agreement with the size distribution obtained from DLS data (black line). Indeed, the DLS data indicate that the peak diameter is 275 ± 120 nm (Figure S6) while the one calculated from electrochemical data is 335 ± 160 nm (Figure S5). The higher value obtained from electrochemical data can be due to the limited sensitivity of the collisions experiments regarding the background current noise which not allow to count the smallest current spikes in the i-t curve. Hence, the calculated liposomes diameter is probably overestimated in the hypothesis where the smallest liposomes (below than 200 nm of diameter) are not detected. In addition, the DLS data indicate that the redox DMPC liposomes solution is quite polydisperse (PDI = 0.19), leading to a quite large size distribution (Figure S6). The results presented in Figure 3 confirm that the current spikes observed in the orange i-t curve in Figure 2 are related to redox liposomes single impacts onto Pt UME. The destructive interaction of RL as a biosurfactant with the liposomes DMPC bilayer is necessary for leading upon impact onto Pt UME to the breakdown of the liposomes and/or the release/electrolysis of its encapsulated redox content. In conclusion, a first proof of concept based on the electrochemistry of redox liposome single impacts for the bacterial virulence factors sensing is proposed with the detection of the lowest concentration of Rhamnolipid (500 nM) in less than 30 minutes. The next step of this study is to optimize the experimental setup such as the electrode and liposomes sizes according to the Rhamnolipid concentration in order to improve the sensitivity of the electrochemical signal and do a statistical analysis relating the time of inoculation to the toxin and liposome concentrations. This work should be extended to other toxins and types of liposomes in order to reach a selective detection which is a crucial parameter for biosensor applications. In the near future, this electrochemical biosensor concept should be used for the detection of virulence factors from the supernatant of pathogenic bacteria.

Entry for the Table of Contents

Insert graphic for Table of Contents here.

A proof of concept showing that the electrochemistry of redox liposomes single impacts is a powerful and efficient analytical technique for sensing target toxins produced by pathogenic bacteria such as Rhamnolipid from Pseudomonas aeruginosa which interacts destructively with the lipid membrane of liposomes. A high sensitivity of sensing has been reached with a submicromolar concentration of Rhamnolipid detected in less than 30 minutes.

Institute and/or researcher Twitter usernames: @CeisamLab

Figure 1 .

 1 Figure 1. Principle of the bacterial Rhamnolipid toxin detection by redox liposome single impact electrochemistry

Figure 2 .

 2 Figure 2. Chronoamperometry measurements recorded at +0.7 V vs Ag/AgCl on 10 µm Pt UME in A. 0.1 M PB at pH 7.4 (black line, left) and B. 500 nM RL and 0.1 M PB at pH 7.4 (black line, right) after addition of 3 pM redox DMPC liposomes at 0 min (blue line), 15 min (orange line), 30 min (green line) and 50 min (pink line).

Figure 3 .

 3 Figure 3. Size distributions from DLS data (black line) and from charge data by integrating current spikes of i-t curve recorded at +0.7 V vs Ag/AgCl on 10 µm Pt UME in 500 nM RL 0.1 M PB at pH 7.4 at 15 min after addition of 3 pM redox DMPC liposomes (orange bar).

Acknowledgements

The support of the Region Pays de la Loire, University of Nantes and French Institute of Sweden are acknowledged. Prof E. Ishow is thank for her helpful training and advices on the DLS instrument. Prof. A. G. Ewing, Prof. J.-M. Bouler and Dr. B. Bujoli are gratefully acknowledged for support and helpful discussions.


Renato Vizuete
email: renato.vizuete@l2s.centralesupelec.fr

Paolo Frasca
email: paolo.frasca@gipsa-lab.fr

Elena Panteley
email: elena.panteley@l2s.centralesupelec.fr

Gradient Descent for Resource Allocation with Packet Loss ⋆

Keywords: Distributed optimization, multi-agent systems, packet loss, gradient descent, resource allocation, open multi-agent systems

This paper studies the effect of packet loss during the application of the weighted gradient descent to solve a resource allocation problem with piecewise quadratic cost functions in a multi-agent system. We define two performance metrics that measure, respectively, the deviation from the constraint and the error on the expected cost function. We derive upper bounds on both metrics: both bounds are proportional to the difference between the initial cost function and the cost function evaluated at the minimizer. Then, we extend the analysis of the constraint violation to open multi-agent systems where agents are replaced: based on a preliminary result and simulations we show that the combination of replacements and losses makes the constraint violation error diverge with time.

INTRODUCTION

The optimal distribution of a budget among multiple entities or activities is one of the most important problems in optimization [START_REF] Ibaraki | Resource Allocation Problems: Algorithmic Approaches[END_REF]. In the case of multi-agent systems, the objective is to find an optimal distribution of the budget among the agents while each agent aims to minimize a local cost function and communications are restricted by the network topology [START_REF] Nedić | Improved convergence rates for distributed resource allocation[END_REF][START_REF] Doan | Distributed resource allocation over dynamic networks with uncertainty[END_REF]. Applications include energy resources [START_REF] Dominguez-Garcia | Decentralized optimal dispatch of distributed energy resources[END_REF], games [START_REF] Liang | Distributed Nash equilibrium seeking for aggregative games with coupled constraints[END_REF], and distributed computer systems [START_REF] Kurose | A microeconomic approach to optimal resource allocation in distributed computer systems[END_REF]. In general, algorithms assume perfect communications and rely on a symmetric exchange of information through the network to maintain the budget fixed (keeping the budget fixed is indeed the constraint of the optimization problem). One of the most important algorithms used to solve the resource allocation problem in a multi-agent system is the weighted gradient descent proposed in [START_REF] Xiao | Optimal scaling of a gradient method for distributed resource allocation[END_REF], where the update rule makes use of a matrix associated to the network that preserves the constraint. However, communication networks can suffer from errors in the communications that cause the loss of information. In this case, even if the graph associated with the communication network was originally undirected, interactions between the agents can effectively be asymmetric when the packet loss occurs only in one direction. The symmetry of the system being broken, relevant quantities may fail to be preserved in many scenarios [START_REF] Patterson | Distributed average consensus with stochastic communication failures[END_REF]Fag-⋆ Research supported in part by the Agence Nationale de la Recherche (ANR) via grant "Hybrid And Networked Dynamical sYstems" (HANDY), number ANR-18-CE40-0010. nani and Zampieri, 2009;[START_REF] Frasca | On the mean square error of randomized averaging algorithms[END_REF]). In the case of the average estimation under packet losses, (Frasca and Hendrickx, 2013a) derived upper bounds for the mean square error to evaluate the deviation from the initial average for general graph topologies. While a possible deviation of the initial average value is not necessarily problematic in the average estimation where the main objective is to reach consensus in most cases, for optimization problems, similar errors can be more consequential, since they may imply the violation of the constraint. For example in the case of the resource allocation problem, agents could perform unnecessary tasks or the sum of the individual tasks may not satisfy the global demand. Indeed, the analysis of the constraint violation in complex scenarios has been the object of many works, specially in the framework of online optimization, where it is assumed that cost functions change in time [START_REF] Yi | Regret and cumulative constraint violation analysis for online convex optimization with long term constraints[END_REF].

Moreover, in large systems, the set of agents can change due to arrivals, departures or replacements, giving rise to open multi-agent systems. Since the size of the system may change with time, it is important to define performance metrics independent of the dimension like the variance (Monnoyer de [START_REF] Monnoyer De Galland | Modelling gossip interactions in open multiagent systems[END_REF] or the normalized mean square error [START_REF] Vizuete | On the influence of noise in randomized consensus algorithms[END_REF]. In the case of optimization in open systems, it is assumed that agents can be replaced during the process, implying a change of the local cost functions [START_REF] Hendrickx | Stability of decentralized gradient descent in open multi-agent systems[END_REF][START_REF] Monnoyer De Galland | Random coordinate descent algorithm for open multi-agent systems with complete topology and homogeneous agents[END_REF]. Due to the possible replacements, the minimizer is time-varying and in many scenarios the original constraint may not be preserved even if communications do not suffer of potential losses.

In this paper, we analyze the effects of packet losses in the performance of the weighted gradient descent algorithm to solve the resource allocation problem. First, assuming possible packet losses, we find conditions to guarantee convergence of the algorithm to the minimizer at least in expectation. Then, we derive upper bounds on the deviation from the constraint and on the error of the expected value of the cost function. Finally, we present an analysis in the case of open multi-agent systems where replacements of the agents can take place and we show, by a preliminary result and simulations, that the constraint error can diverge for this particular scenario. This fact highlights the difference between closed and open multiagent systems with packet loss.

PROBLEM FORMULATION

The set of real numbers is denoted by R and the set of non negative integers by Z ≥0 . For two vectors x, y ∈ R n , we denote the Euclidean inner product by ⟨x, y⟩ = x T y and the Euclidean norm by ||x|| = (x T x) 1/2 . The vector constituted of only ones is denoted by 1. The Hadamard or element-wise product between two matrices A, B ∈ R n×n is denoted by A ⊙ B.

Resource Allocation Problem

The resource allocation problem in a system of n agents is formulated as:

min x=[x1,...,xn]∈R n f (x) = n i=1 f i (x i ) subject to n i=1 x i = b, (1)
where f i : R → R is the local cost function associated to agent i and b is the budget. Following standard assumptions, we consider that each local cost function is continuously differentiable, β-smooth (i.e., |f ′ i (x) -f ′ i (y)| ≤ β |x -y| , ∀x, y) and α-strongly convex (i.e., f i (x) -α 2 |x| 2 is convex)(see [START_REF] Nesterov | Lectures on Convex Optimization[END_REF]). Hence, the global cost function f : R n → R is also β-smooth and α-strongly convex and the optimal solution denoted by x * is unique. Furthermore, the optimality condition implies that 1 T x * = b and ∇f * = ζ1 for some ζ ∈ R where f * = f (x * ).

We consider that agents interact through a connected undirected graph G = (V, E), where V = {1, . . . , n} is the set of agents and E ⊆ V × V is the set of edges. To the graph G we associate an adjacency matrix A = [a ij ] where a ij = 1 if (i, j) ∈ E and a ij = 0 otherwise. We assume that the graph does not have self-loops (i.e., a ii = 0). The Laplacian matrix of the graph G is defined as L = D -A where D = diag[d 1 , . . . , d n ] and d i is the degree of agent i (i.e., d i is the ith row-sum of A). Notice that L is a symmetric matrix that satisfies L = L T .

To solve problem (1) at each iteration we apply the weighted gradient descent algorithm with positive stepsize h:

x t+1 = x t -hL∇f (x t ),
(2) where L determines the exchange of information between the agents according to the network topology. It is wellknown that for any initial condition x 0 satisfying the constraint 1 T x 0 = b, algorithm (2) solves problem (1) with an appropriate choice of the step-size h and guarantees linear convergence of the cost function.

Packet Losses

In this paper, we consider that the communication network is unreliable and some messages can be lost. Following the approach in (Frasca and Hendrickx, 2013a) we assume that at each time instant t ∈ Z ≥0 , we have a matrix

Γ t = γ ij t ∈ {0, 1} n×n composed of Bernoulli random
variables γ ij t independent across t, i, j with P(γ ij t = 1) = p for some p ∈ (0, 1). Under this assumption, a link a ij = 1 is active if γ ij t = 1 and inactive if γ ij t = 0. If a ij = 0, the value of γ ij t does not have any effect. The topology of the network with possible packet losses at each time instant is thus defined by:

A t = A ⊙ Γ t , (3)
L t = D t -A t ,
(4) and the weighted gradient descent algorithm becomes

x t+1 = x t -hL t ∇f (x t ).
(5) Due to the asymmetric loss of information, the graph topology at each time instant is directed and hence the Laplacian matrix L t is non-symmetric. The constraint can be violated during an update of the states of the network using ( 5) and hence, the minimizer of the problem (1) will not be reached.

At least, the algorithm (5) must solve (1) in expectation. Since Γ t is independent in time and E [L t ] = pL, the expected dynamics is given by:

E [x t+1 ] = E [x t ] -hpLE [∇f (x t )] .
(6) Notice that similar to (2), the dynamics (6) preserves the constraint in expectation since:

1

T E [x t+1 ] = 1 T (E [x t ] -hpLE [∇f (x t )]) = 1 T E [x t ] . (
7) However, the minimizer x * is not necessarily a stationary point of (6). Indeed, assume that E [x t ] = x * , then we have

E [x t+1 ] = x * -hLE [∇f (x t )] .
Under the assumption that the local functions are βsmooth and α-strongly convex we cannot guarantee that

E [∇f (x t )] = ∇ (f (E [x t ])) = ζ1.
For this reason, we restrict the set of possible cost functions to piecewise quadratic functions, which are β-smooth and α-strongly convex. Assumption 1. Each local cost function f i (x i ) is continuous differentiable piecewise quadratic (i.e., there exists finitely many quadratic functions {q j } J j=1 such that

f i (x i ) ∈ {q j (x i )} J j=1 for all x i ∈ R).
Remark that under Assumption 1, the gradient is a vector of piecewise linear functions [∇f (x)] i = ∂f (x) ∂xi = m ij x i +c ij and we have

E [∇f (x t )] = ∇f (E [x t ]) .
(8) Then, the expected dynamics (6) becomes:

E [x t+1 ] = E [x t ] -hpL∇f (E [x t ]) ,
(9) which corresponds to the original weighted gradient descent, thereby guaranteeing that the constraint is preserved at each time instant and that the minimizer is reached such that lim

t→∞ E [x t ] = x * and lim t→∞ f (E [x t ]) = f * .
However, a similar equality is not true for the expected value of the cost function. Since the functions are convex, by using Jensen's inequality we obtain:

E [f (x t )] ≥ f (E [x t ]
) , and by taking the limit (whose existence will be discussed in the next section) we get:

lim t→∞ E [f (x t )] ≥ f * , ( 10
)
where the inequality is in general strict.

Performance of the gradient descent algorithm (5) in the presence of packet losses can be characterized by the following two quantities: errors due to the violation of the constraint and errors due to the deviation of the expected cost from the ideal cost f * . In the following sections, we analyze these two errors by defining appropriate performance metrics.

CONSTRAINT VIOLATION

Although the constraint is preserved in expectation according to (7), this fact does not imply that the constraint is also preserved for any particular realization of the process. In this section we estimate the deviation from the constraint by deriving an upper bound for the constraint violation metric:

J constr := lim sup t→∞ E 1 T x t -b 2 , (11)
where the initial condition x 0 satisfies 1 T x 0 = b. Notice that (11) corresponds to the asymptotic mean square error.

Theorem 1. Consider the weighted gradient descent given by ( 3), ( 4) and ( 5). Under Assumption 1, for any positive scalar h ≤ 2 β(pλn+2-2p) , the constraint violation metric satisfies:

J constr ≤ 4h(1 -p) 2 -hβ(pλ n -2p + 2) (f (x 0 ) -f * ) , ( 12
)
where λ n is the largest eigenvalue of L.

Before presenting the proof of Theorem 1, we recall a lemma derived in (Frasca and Hendrickx, 2013a). Lemma 1. For L t given by ( 3) and ( 4):

E L T t L t = p 2 L 2 + 2p(1 -p)L (13) E L T t 11 T L t = 2p(1 -p)L. ( 14
) Proof of Theorem 1. Let us denote H t = 1 T x t -b .
For (5) we denote by {F t } t∈Z ≥0 the filtration of σ-algebras generated by the process x t . We compute the expectation conditioned upon the filtration generated by x t :

E H 2 t+1 |F t = E 1 T (x t -hL t ∇f (x t )) -b 2 = H 2 t + h 2 ∇f (x t ) T E L T t 11 T L t ∇f (x t ) -hH t 1 T E [L t ] ∇f (x t ) = H 2 t + h 2 ∇f (x t ) T E L T t 11 T L t ∇f (x t )
, where we have used the fact that 1 T L = 0. Then, we use (14) to obtain:

E H 2 t+1 |F t = H 2 t + 2h 2 p(1 -p) L 1/2 ∇f (x t ) 2 .
If we take the total expectation we get:

E E H 2 t+1 |F t =E H 2 t +2h 2 p(1-p)E L 1/2 ∇f (x t ) 2 .
By using a recursive argument we have that at the time instant t:

E H 2 t = H 2 0 + 2h 2 p(1 -p) t-1 s=0 E L 1/2 ∇f (x s ) 2 = 2h 2 p(1 -p) t-1 s=0 E L 1/2 ∇f (x s ) 2 , ( 15
)
where we used the fact that at the initial time t = 0, the constraint is not violated H 0 = 0. Now, let us consider the inequality corresponding to β-smooth functions:

f (x t+1 ) ≤ f (x t ) + β 2 ||x t -x t+1 || 2 + ⟨∇f (x t ), x t+1 -x t ⟩ = f (x t ) + h 2 β 2 ⟨∇f (x t ), L T t L t ∇f (x t )⟩ -h⟨∇f (x t ), L t ∇f (x t )⟩.
We compute the expectation given the filtration generated by x t and we use (13) to get:

E [f (x t+1 )|F t ] ≤ f (x t ) -h⟨∇f (x t ), pL∇f (x t )⟩ + h 2 β 2 ⟨∇f (x t ), (p 2 L 2 + 2p(1 -p)L)∇f (x t )⟩.
We use the upper bound ||L∇f

(x t )|| 2 ≤ λ n L 1/2 ∇f (x t )
2 and we obtain:

E [f (x t+1 )|F t ] ≤ f (x t ) -hp 1 -hβ(1 -p) - hβpλ n 2 L 1/2 ∇f (x t ) 2 .
We denote ξ = hp 1 -hβ(1 -p) -hβpλn 2 and we compute the total expectation to get:

E [E [f (x t+1 )|F t ]] ≤ E [f (x t )] -ξE L 1/2 ∇f (x t ) 2 . (16) Since h ≤ 2 β(pλn+2-2p) we guarantee that E [f (x t+1 )] ≤ E [f (x t )] such that: E L 1/2 ∇f (x t ) 2 ≤ 1 ξ (E [f (x t )] -E [f (x t+1 )]) .
If we take the sum over the time instants s we have:

t-1 s=0 E L 1/2 ∇f (x s ) 2 ≤ 1 ξ (f (x 0 ) -E [f (x t )]) . Notice that since E [f (x t )] is monotonically decreasing and E [f (x t )] ≥ f (x * )
, where x * is the minimizer of f (x) without constraints, the limit lim t→∞ E [f (x t )] exists and we get:

lim t→∞ t-1 s=0 E L 1/2 ∇f (x s ) 2 ≤ 1 ξ f (x 0 ) -lim t→∞ E [f (x t )] ≤ 1 ξ (f (x 0 ) -f * ) , (17)
where we used the fact that

E [f (x t )] ≥ f (E [x t ]
). Finally, by using ( 17) in (15), we complete the proof. 2

By inspecting the right-hand side of (12), we can make several observations. The upper bound decreases when p increases and lim

p→1 J constr = 0, (18)
which is consistent with the fact that the constraint is always preserved with no communication losses. In the case of multiple failures, which correspond to p small, we parametrize the step size as h = 2ϵ β(pλn+2-2p) with ϵ ∈ (0, 1) and using ( 12) we obtain:

lim p→0 J constr ≤ 2ϵ β(1 -ϵ) (f (x 0 ) -f * ) . (19)
From ( 19), we see that for small values of p, large values of ϵ, corresponding to large values of the step-size, increase the potential violation of the constraint. Figure 1 presents the computation of the upper bound for the constraint violation metric in a complete graph with n = 10 agents, β = 1 and ϵ = 0.5, where we can observe the behavior of the upper bound and the limit values determined by ( 18) and ( 19).

ERROR ON THE EXPECTED COST FUNCTION

Due to the packet losses, f (x t ) is a random process whose evolution depends on the different realizations of L t . When the degraded version of the weighted gradient descent algorithm ( 5) is executed, it is important to know the distribution of f (x t ), to evaluate the impact of the perturbations on the performance of ( 5) and identify the additional cost corresponding to the deviation from the ideal cost f * . From ( 10) we can observe that E [f (x t )] is greater than f * . Since the objective is to remain as close as possible to f * , we define the cost function metric:

J funct := lim sup t→∞ E [f (x t )] -f (x * ), ( 20
)
where the initial condition x 0 satisfies 1 T x 0 = b. The metric (20) estimates the gap in the cost function due to the packet losses during the application of the weighted gradient descent algorithm. It is important to notice that (20) corresponds to the limit of the so called Jensen gap [START_REF] Abramovich | Some new estimates of the 'Jensen gap[END_REF].

Theorem 2. Consider the weighted gradient descent given by ( 3), ( 4) and ( 5). Under Assumption 1, for any positive scalar h ≤ 2 β(pλn+2-2p) , the cost function metric satisfies:

J funct ≤ 1 - 2 -hβ (pλ n -2p + 2) 2 -hαpλ 2 (f (x 0 ) -f * ) , (21
) where λ 2 is the second smallest eigenvalue of L.

Proof. From ( 16), we have that the total expectation of f (x t+1 ) satisfies:

E [f (x t+1 )] ≤ E [f (x t )] -ξE L 1/2 ∇f (x t ) 2 ,
and by using Jensen's inequality and (8) we get:

E [f (x t+1 )] ≤ E [f (x t )] -ξ E L 1/2 ∇f (x t ) 2 = E [f (x t )] -ξ L 1/2 ∇f (E [x t ]) 2 .
If we take the sum over the time instants s we obtain:

E [f (x t )] ≤ f (x 0 ) -ξ t-1 s=0 L 1/2 ∇f (E [x s ]) 2 ,
and the limit satisfies:

lim t→∞ E [f (x t )] ≤ f (x 0 ) -ξ lim t→∞ t-1 s=0 L 1/2 ∇f (E [x s ]) 2 .
(22) Now, we consider the dynamics in expectation given by ( 6), which satisfies the inequality of α-strongly convex functions:

f (E [x t+1 ]) ≥ f (E [x t ]) + ⟨∇f (E [x t ]), E [x t+1 ] -E [x t ]⟩ + α 2 ||E [x t+1 ] -E [x t ]|| 2 = f (E [x t ]) -hp L 1/2 ∇f (E [x t ]) 2 + h 2 αp 2 2 ||L∇f (E [x t ])|| 2 .
Notice that due to the Min-Max Theorem, for any vector z ∈ R n , it holds

||Lz|| 2 = L 1/2 L 1/2 z 2 ≥ λ 2 L 1/2 z 2 , (23)
since L 1/2 z is orthogonal to 1. Hence we obtain:

f (E [x t+1 ]) ≥ f (E [x t ]) -hp 1-hαpλ2 2 L 1/2 ∇f (E [x t ]) 2 . We denote δ = hp 1 -hαpλ2 2 , which is non-negative if h ≤ 2 αλ2 . Since 2 β(pλn+2-2p) ≤ 2 αλ2 we get: L 1/2 ∇f (E [x t ]) 2 ≥ 1 δ (f (E [x t ]) -f (E [x t+1 ])) .
If we take the sum over the time instants s we obtain:

t-1 s=0 L 1/2 ∇f (E [x s ]) 2 ≥ 1 δ (f (x 0 ) -f (E [x t ])) ,
and for p > 0 the limit satisfies:

lim t→∞ t-1 s=0 L 1/2 ∇f (E [x s ]) 2 ≥ 1 δ (f (x 0 ) -f * ) . (24)
We use ( 24) in ( 22) to get:

lim t→∞ E [f (x t )] ≤ f (x 0 ) - ξ δ (f (x 0 ) -f * ) , which yields lim t→∞ E [f (x t )] -f (x * ) ≤ 1 - ξ δ (f (x 0 ) -f * ) . 2
To analyze the asymptotic behavior of the upper bound for J funct , let us consider a step-size given by h = 2ϵ β(pλn+2-2p) with ϵ ∈ (0, 1). For low values of probabilities corresponding to non-robust networks we have: lim which implies that small step-sizes, characterized by small values of ϵ, reduce more the gap between the expected value of the cost function and f * . In the case of large values of probabilities, we have the following asymptotic behavior:

p→0 J funct ≤ ϵ (f (x 0 ) -f * ) , (25)
lim p→1 J funct ≤ ϵ(1 -κ -1 κ -1 λ ) 1 -ϵκ -1 κ -1 λ (f (x 0 ) -f * ) , (26)
where κ = β α is the condition number of the cost functions f i and κ λ = λn λ2 . We can observe that only for the particular case of a complete graph with κ = 1 we obtain lim

p→1 J funct = 0, (27)
which implies a certain level of conservatism in the upper bound since ( 27) should hold for any choice of κ and κ λ in order to recover the behavior of a network without packet losses. Figure 2 shows the computation of the upper bound for J funct for a complete graph and a star graph with n = 10 agents, ϵ = 0.5, f (x 0 ) -f * = 1 and two values of κ where we can appreciate the behavior of the upper bound and the limit values defined by ( 25) and (26).

OPEN MULTI-AGENT SYSTEMS

Let us consider the open scenario when agents can be replaced at any time step. In order to preserve the constraint of problem (1) at least in expectation, we shall we consider only replacements such that the local cost function of an agent f i can change but the state x i is kept. More precisely, we assume that at each time instant t we can have an update with probability p U or a replacement with probability p R = 1 -p U . To simplify notations, we will denote f t , the global cost function at the time instant t. If an update occurs, the weighted gradient descent is applied and f t+1 = f t . If a replacement occurs, a single agent i is uniformly randomly selected and receives a new (piecewise quadratic) cost function.

When the open system does not suffer from packet losses, the constraint is always preserved but the global minimizer x * t is time-varying since the cost functions change in time. In this analysis, we assume that the open system also suffers of packet losses such that the weighted gradient descent algorithm becomes:

x t+1 = x t -hL t ∇f t (x t ).

(28) Similarly to the closed system we analyze the constraint violation metric J constr defined in (11).

Proposition 1. Consider the weighted gradient descent given by ( 3), ( 4) and ( 28). Under Assumption 1, for any positive scalar h ≤ 2 β(pλn+2-2p) , the constraint violation metric satisfies:

J constr ≥ 4p U h(1-p) 2-αh(pλ2+2-2p) lim t→∞ t-1 s=0 E [f s (x s ) -f s (x s+1 )] . (29)
Proof. As in the proof of Theorem 1, we use the notation for the error H t = 1 ⊤ x t -b and the natural filtration {F t } t∈Z ≥0 . If at the moment t, an update U happens, then from ( 15) we have at t + 1:

E H 2 t+1 |U = E H 2 t + 2h 2 p(1 -p)E L 1/2 ∇f t (x t ) 2 .
(30) When a replacement R occurs, the constraint is not violated at that event

E H 2 t+1 |R = E H 2 t .
Then, the total expectation at the time instant t + 1 is defined as:

E H 2 t+1 = p U E H 2 t+1 |U + p R E H 2 t+1 |R = E H 2 t + 2p U h 2 p(1 -p)E L 1/2 ∇f t (x t ) 2 .
Using a recursive argument we obtain:

E H 2 t = 2p U h 2 p(1 -p) t-1 s=0 E L 1/2 ∇f s (x s ) 2 . (31)
Since the cost functions f t are α-strongly convex we have:

f t (x t+1 ) ≥ f t (x t ) + α 2 ||x t -x t+1 || 2 + ⟨∇f t (x t ), x t+1 -x t ⟩ = f t (x t ) + h 2 α 2 ⟨∇f t (x t ), L T t L t ∇f t (x t )⟩ -h⟨∇f t (x t ), L t ∇f t (x t )⟩.
We compute the expectation given the filtration generated by x t and we use (13) to get:

E [f t (x t+1 )|F t ] ≥ f t (x t ) -h⟨∇f t (x t ), pL∇f t (x t )⟩ + h 2 α 2 ⟨∇f t (x t ), (p 2 L 2 + 2p(1 -p)L)∇f t (x t )⟩.
By using (23), we have the lower bound ||L∇f

t (x t )|| 2 ≥ λ 2 L 1/2 ∇f t (x t )
2 and we obtain:

E [f t (x t+1 )|F t ] ≥ f t (x t ) -hp 1 -hα(1 -p) - hαpλ 2 2 L 1/2 ∇f t (x t ) 2 . We denote θ = hp 1 -hα(1 -p) -hαpλ2 2
and we compute the total expectation to get:

E [E [f t (x t+1 )|F t ]] ≥ E [f t (x t )] -θE L 1/2 ∇f t (x t ) 2 . Since h ≤ 2 β(pλn+2-2p) ≤ 2 α(pλ2+2-2p) we guarantee that E [f t (x t+1 )] ≤ E [f t (x t )] and therefore E L 1/2 ∇f t (x t ) 2 ≥ 1 θ (E [f t (x t )] -E [f t (x t+1 )]) . ( 32
)
By using (32) in (31) we complete the proof. 2 Unlike a closed system, we observe that due to the possible replacements that take place in the system, in general . This mismatch implies that terms are not canceled in the series in the lower bound and the error accumulates in time, so that J constr grows to infinity. This intuition is corroborated by the simulations of E (1 T x t -b) 2 for a complete graph with n = 7 agents that are shown in Figure 3. We can observe in the plot that both the constraint violation metric and the lower bound diverge. When p = 1, corresponding to an ideal system without packet losses, even if replacements happen in the system (i.e., p U > 0), J constr is zero since E H 2 t+1 |U = E H 2 t as we can see in (30). This can also be derived from ( 7), since the algorithm preserves the constraint independently of the cost functions.

CONCLUSION

In this paper, we analyzed the effect of packet losses in the performance of the weighted gradient descent algorithm. We defined two performance metrics to measure the deviation from the constraint and the error of the expected cost function and we derived an upper bound. We extended the analysis of packet losses to open multi-agent systems and we showed that the constraint violation metric may diverge in such a case.

Several questions remain open as future work. First, it would be useful to improve the upper bound for the cost function metric since the current result suffers from a certain level of conservatism. Second, the analysis of open case is very preliminary, even though it has already shown that its properties are very different from those of the closed case. Relevant questions in this direction include the improvement of the lower bound for the constraint violation metric and the analysis of the cost function metric. Finally, it would be interesting to extend the analysis of this work to algorithms based on the computation of the gradient of a subset of agents like the random coordinate descent algorithm [START_REF] Necoara | Random coordinate descent algorithms for multi-agent convex optimization over networks[END_REF][START_REF] Monnoyer De Galland | Random coordinate descent algorithm for open multi-agent systems with complete topology and homogeneous agents[END_REF].

Fig. 1 .

 1 Fig. 1. Computation of the upper bound (12) for J constr for a complete graph with n = 10, β = 1, ϵ = 0.5 and f (x 0 ) -f * = 1.

Fig. 2 .

 2 Fig. 2. Computations of the upper bound (21) for J funct for a complete and star graph with n = 10, ϵ = 0.5, f (x 0 ) -f * = 1 and two values of κ.

Fig. 3 .

 3 Fig. 3. Simulation of E (1 T x t -b) 2 and of the lower bound (29) for a complete graph under replacements with n = 7, α = 1, β = 10, p = 0.8, p U = 0.5 by considering 10000 realizations of the process.E [f s (x s+1 )] ̸ = E [f s+1 (x s+1 )]. This mismatch implies that terms are not canceled in the series in the lower bound and the error accumulates in time, so that J constr grows to infinity. This intuition is corroborated by the simulations of E (1 T x t -b) 2 for a complete graph with n = 7 agents that are shown in Figure3. We can observe in the plot that both the constraint violation metric and the lower bound diverge. When p = 1, corresponding to an ideal system without packet losses, even if replacements happen in the system (i.e., p U > 0), J constr is zero sinceE H 2 t+1 |U = E H 2