Alexis Bordet

Ryan F Landis

Yiwei Lee

Juan Manuel Tonga

Chenghsuan Asensio

Pier-Francesco Li

Aikaterini Fazzini

Vincent M Soulantika

Bruno Rotello
email: rotello@chem.umass.edu

Chaudret
email: chaudret@insa-toulouse.fr

Gulen Y Tonga

Juan M Asensio

Chenghsuan Li

Pier-Francesco Fazzini

Katerina Soulantica

Vincent M Rotello

Bruno Chaudret

Water-Dispersible and Biocompatible Iron Carbide Nanoparticles with High Specific Absorption Rate

published or not. The documents may come   L'archive ouverte pluridisciplinaire

cellular uptake of iron carbide nanoparticles could be easily tuned and were highly dependent on the chemical structure of the ligands used. KEY WORDS. nanoparticles -iron carbide -magnetic hyperthermia -functionalizationbiocompatibility Magnetic nanoparticles (NPs) are of strategic interest for several biomedical applications such as magnetic resonance imaging, 1-3 magnetic hyperthermia [START_REF] Fernandes Cardoso | Advances in Magnetic Nanoparticles for Biomedical Applications[END_REF][START_REF] Lee | Exchange-Coupled Magnetic Nanoparticles for Efficient Heat Induction[END_REF][START_REF] Salunkhe | Magnetic Hyperthermia with Magnetic Nanoparticles: a Status Review[END_REF][START_REF] Kumar | Magnetic Nanomaterials for Hyperthermia-Based Therapy and Controlled Drug Delivery[END_REF] and drug delivery. [START_REF] Fernandes Cardoso | Advances in Magnetic Nanoparticles for Biomedical Applications[END_REF][START_REF] Sun | Magnetic Manoparticles in MR Imaging and Drug Delivery[END_REF][START_REF] Kumar | Magnetic Nanomaterials for Hyperthermia-Based Therapy and Controlled Drug Delivery[END_REF] When exposed to an alternating magnetic field, suitable magnetic NPs undergo energy losses, thus producing heat. This phenomenon, called magnetic hyperthermia, is widely applied biomedically to destroy or weaken tumors, [START_REF] Salunkhe | Magnetic Hyperthermia with Magnetic Nanoparticles: a Status Review[END_REF] trigger drug delivery [START_REF] Kumar | Magnetic Nanomaterials for Hyperthermia-Based Therapy and Controlled Drug Delivery[END_REF] or improve drug efficacy. [START_REF] Kumar | Magnetic Nanomaterials for Hyperthermia-Based Therapy and Controlled Drug Delivery[END_REF] For magnetic hyperthermia to be efficient, however, magnetic NPs must be designed to maximize energy losses and thus the production of heat. The ability of NPs to heat through magnetic induction is commonly expressed by their specific absorption rate (SAR), and the objective is to reach high SAR at frequencies and amplitudes harmless for healthy tissues [START_REF] Giustini | Magnetic Nanoparticle Hyperthermia in Cancer Treatment[END_REF][START_REF] Hergt | Magnetic Particle Hyperthermia-Biophysical Limitations of a Visionary Tumour Therapy[END_REF] (typically 100 kHz, 20 mT). Substantial effort has been dedicated to the synthesis of water stable and biocompatible NPs displaying high SARs. Given these requirements, iron-based nanoparticles appear to be promising candidates. [START_REF] Hergt | Magnetic Particle Hyperthermia: Nanoparticle Magnetism and Materials Development for Cancer Therapy[END_REF] Among them, iron oxide nanoparticles (IONPs) have been extensively studied. As a result, a significant number of synthetic procedures allowing the production of tunable iron oxide nanoparticles (size, shape, composition, etc.) are available. [START_REF] Gupta | Synthesis and Surface Engineering of Iron Oxide Nanoparticles for Biomedical Applications[END_REF][START_REF] Laurent | Magnetic Iron Oxide Nanoparticles: Synthesis, Stabilization, Vectorization, Physicochemical Characterizations, and Biological Applications[END_REF] While commercially available IONPs generally display poor hyperthermia properties, very good SARs were reported for complex IONPs such as magnetite bacterial magnetosomes, [START_REF] Hergt | Magnetic Properties of Bacterial Magnetosomes As Potential Diagnostic and Therapeutic Tools[END_REF] maghemite nanoflowers [START_REF] Hugounenq | Iron Oxide Monocrystalline Nanoflowers for Highly Efficient Magnetic Hyperthermia[END_REF] or maghemite nanocubes. [START_REF] Guardia | Water-Soluble Iron Oxide Nanocubes with High Values of Specific Absorption Rate for Cancer Cell Hyperthermia Treatment[END_REF][START_REF] Guardia | One Pot Synthesis of Monodisperse Water Soluble Iron Oxide Nanocrystals with High Values of the Specific Absorption Rate[END_REF] Despite their promising properties, these complex IONPs are not easily synthesized and exhibit a low saturation magnetization (70-80 A.m 2 .kg -1 in the best cases). Iron has the highest saturation magnetization (210 A.m 2 .kg -1 ) among all elements, and Fe(0) NPs were shown to have very good SARs in organic solvents. [START_REF] Mehdaoui | Large Specific Absorption Rates in the Magnetic Hyperthermia Properties of Metallic Iron Nanocubes[END_REF][START_REF] Mehdaoui | Large Specific Absorption Rates in the Magnetic Hyperthermia Properties of Metallic Iron Nanocubes[END_REF] Unfortunately, in the presence of oxygen or water, zero-valent iron NPs undergo very fast oxidation that dramatically affects their magnetic properties, reducing their utility as hyperthermic agents.

Iron carbide NPs provide a potential compromise between high saturation magnetization and stability in biological media. Several studies evidenced their potential for applications including magnetic hyperthermia, [START_REF] Yu | Iron Carbide Nanoparticles: an Innovative Nanoplatform for Biomedical Applications[END_REF][START_REF] Meffre | A Simple Chemical Route toward Monodisperse Iron Carbide Nanoparticles Displaying Tunable Magnetic and Unprecedented Hyperthermia Properties[END_REF][START_REF] Bordet | Magnetically Induced Continuous CO2 Hydrogenation Using Composite Iron Carbide Nanoparticles of Exceptionally High Heating Power[END_REF] magnetic resonance imaging, [START_REF] Tang | Fe5C2 Nanoparticles with High MRI Contrast Enhancement for Tumor Imaging[END_REF] catalysis. [START_REF] Meffre | A Simple Chemical Route toward Monodisperse Iron Carbide Nanoparticles Displaying Tunable Magnetic and Unprecedented Hyperthermia Properties[END_REF][START_REF] Bordet | Magnetically Induced Continuous CO2 Hydrogenation Using Composite Iron Carbide Nanoparticles of Exceptionally High Heating Power[END_REF][START_REF] Yang | Fe5C2 Nanoparticles: A Facile Bromide-Induced Synthesis and As an Active Phase for Fischer-Tropsch Synthesis[END_REF][START_REF] Niether | Improved Water Electrolysis Using Magnetic Heating of Core-Shell FeC@Ni Nanoparticles[END_REF][START_REF] Bordet | Enhancement of Carbon Oxides Hydrogenation on Iron-Based Nanoparticles by In-Situ Water Removal[END_REF] We recently reported the synthesis of specific Fe2.2C ICNPs displaying exceptional SAR for applications in CO2 hydrogenation. [START_REF] Bordet | Magnetically Induced Continuous CO2 Hydrogenation Using Composite Iron Carbide Nanoparticles of Exceptionally High Heating Power[END_REF] We report here the use of dopamine-based ligands to functionalize these ICNP structures and stabilize them efficiently in aqueous and biological media. Even after extended exposure to water and oxygen, these ICNPs present SARs higher than the ones of Fe(0) nanoparticles. The critical influence of the ligand structural design regarding the cytotoxicity and cellular-uptake of these ICNPs was also evidenced.

RESULTS/DISCUSSION

Monodisperse ICNPs (15 nm) were synthesized through the carbidization of preformed Fe(0) nanoparticles under CO/H2, following an approach previously described by our group (Figure S1). [START_REF] Bordet | Magnetically Induced Continuous CO2 Hydrogenation Using Composite Iron Carbide Nanoparticles of Exceptionally High Heating Power[END_REF][START_REF] Bordet | A New Approach to the Mechanism of Fischer-Tropsch Syntheses Arising from Gas Phase NMR and Mass Spectrometry[END_REF] These nanoparticles displayed excellent hyperthermia properties, with specific absorption rates (SARs) far higher than any iron oxide or iron(0) nanoparticles. However, due to the fatty acid and amine ligands of the as-prepared ICNPs, these particles are dispersible only in organic solvents, and further surface modification was mandatory to provide solubility and colloidal stability in aqueous environment. This challenge was addressed by replacing the initial hydrophobic stabilizers by specifically designed dopamine-based ligands (Figure 1) through a ligand exchange reaction. The ligand family was synthesized by adapting a previously described synthetic procedure [START_REF] Miranda | Enzyme-Amplified Array Sensing of Proteins in Solution and in Biofluids[END_REF][START_REF] Yan | Characterization of Surface Ligands on Functionalized Magnetic Nanoparticles Using Laser Desorption/Ionization Mass Spectrometry (LDI-MS)[END_REF] starting from the dopamine derivative (I), which contains a catechol moiety. This catechol moiety is known to bind strongly to iron-based nanoparticle surfaces, [START_REF] Hong | Surface PEGylation and Ligand Exchange Chemistry of FePt Nanoparticles for Biological Applications[END_REF] and especially to iron oxide. However, dopamine derivatives were also recently described as suitable ligands for the stabilization of iron carbide nanoparticles (Fe5C2) in water. [START_REF] Yu | Iron Carbide Nanoparticles: an Innovative Nanoplatform for Biomedical Applications[END_REF][START_REF] Tang | Fe5C2 Nanoparticles with High MRI Contrast Enhancement for Tumor Imaging[END_REF][START_REF] Yu | Multistimuli-Regulated Photochemothermal Cancer Therapy Remotely Controlled via Fe5C2 Nanoparticles[END_REF] The ligand structure contains a hydrophobic alkyl chain (II) which contributes to the stabilization of the NPs by fostering micelle-like hydrophobic assembly of the ligands. The tetraethylene glycol moiety (III) affords water solubility, biocompatibility and stability to the functionalized NPs. At the terminus, the function (IV) of the ligand can be easily modified to tune the behavior of the NPs in biological media. In a typical ligand exchange reaction, ICNPs (10 mg) were dispersed in a mixture of dichloromethane and methanol (5/1). The water-soluble ligand (50 mg, 5 mass equivalents compared to ICNPs) was added to the dispersion, followed by triethylamine (250 µL). The mixture was stirred overnight at room temperature. The work up of the reaction consisted of 3 washings with hexane followed by redispersion in water and purification by dialysis for 48 h. The presence of the desired ligands on the surface of ICNPs was confirmed by LDI/MS experiments (Figure S2-6). [START_REF] Yan | Characterization of Surface Ligands on Functionalized Magnetic Nanoparticles Using Laser Desorption/Ionization Mass Spectrometry (LDI-MS)[END_REF] In addition, thermo-gravimetric analysis of the ICNPs before and after the ligand exchange reaction was performed (Figure S7).

The results showed that the characteristic signature of the long chain amine and acid ligands observed for as-synthesized ICNPs is missing in the case of water soluble ICNPs, and replaced by a signal attributed to the dopamine-based ligands. These observations indicate that the long chain amine and acid ligands were efficiently exchanged by the dopamine-derivatives during the ligand exchange reaction. This is in agreement with LDI-MS measurements, which did not evidence the presence of the initial hydrophobic ligands at the surface of the functionalized ICNPs (Figure S2-6). Taken together, these results suggest that the ICNPs are stabilized through the direct coordination of the dopamine-based ligands, rather than by their intercalation through hydrophobic interaction with remaining hexadecylamine and palmitic acid.

The stability of ICNPs functionalized with different types of ligands was investigated in multiple aqueous media (de-ionized water, phosphate buffer 5 mM (PB), phosphate-buffered saline 150 mM (PBS) and media containing 10% bovine serum albumin (BSA) using dynamic light scattering (DLS) and transmission electron microscopy (TEM) analysis. Special interest was dedicated to the stability of ICNPs in PBS and serum media since their properties are closer to real biological conditions. The multiple types of ligands tested -commercial or synthesized -are disclosed in Table 1, together with their ability to bring ICNPs solubility and stability in aqueous media. The first stability tests were performed with simple commercial ligands commonly used to stabilize magnetic NPs such as iron oxide in aqueous media. [START_REF] Korpany | Stable Water-Soluble Iron Oxide Nanoparticles Using Tiron[END_REF][START_REF] Liu | Facile Surface Functionalization of Hydrophobic Magnetic Nanoparticles[END_REF][START_REF] Wei | Compact Zwitterion-Coated Iron Oxide Nanoparticles for In Vitro and In Vivo Imaging[END_REF] The results show that these ligands all fail to stabilize ICNPs in water, with the nanoparticles aggregating and precipitating in a few minutes. ICNPs thus require more advanced ligands than iron oxide NPs to be stable in water. The use of specifically synthesized dopamine-based ligands provided ICNPs with surprisingly good solubility and stability in all aqueous media studied. Among the five dopamine-based ligands presented in this table, we will focus on three particular ligands displaying respectively anionic (DOP-TEG-COOH), cationic (DOP-TEG-C6) and zwitterionic functionality (DOP-TEG-Zwitter).

For all three ligands, TEM analyses of the NPs before and after functionalization show that the NPs morphology is not affected by the ligand exchange or water exposure (Figure 2a and2c).

Interestingly, the stability of the functionalized ICNPs is such that they cannot be precipitated by magnet-assisted decantation. At suitable concentration, however, the solutions containing functionalized ICNPs exhibit a ferrofluid behavior (Figure 2d). According to DLS measurements, the ICNPs display very good colloidal stability in all types of aqueous media considered, with their hydrodynamic diameter remaining constant for more than one month (Figure 2e,g,i). While in all cases ICNPs appear to be individually stabilized in water, PB and PBS (DH ~ 15-20 nm), the hydrodynamic diameter of the particles in serum media is significantly larger (~ 90 nm). This increase in size is attributed to a partial aggregation of nanoparticles induced by the interactions between the ligands and the proteins of the serum, forming a protein corona shell. [START_REF] Del Pino | Protein Corona Formation Around Nanoparticles -From the Past To the Future[END_REF] TEM analysis revealed the presence of ~ 90 nm assemblies, thus confirming this hypothesis (Figure S8). XRD analysis evidences partial oxidation of the nanoparticles after a four-month exposure to water and air, but the pseudo-hexagonal Fe2.2C structure is still clearly visible (Figure S9). As a result of oxidation, VSM analysis (Figure 3a) shows that the saturation magnetization (90 A.m 2 .kg -1 ) and the coercive field (32 mT) of ICNPs are significantly reduced as compared to starting nanoparticles (160 A.m 2 .kg -1 , 102 mT respectively). Nevertheless, the SAR of these functionalized NPs in water remain high, reaching ca. 1000 W•g -1 at 47 mT, 100 kHz. Indeed, their magnetic heating properties even after partial oxidation are still better than those of 12.5 nm Fe(0) NPs synthesized and kept under oxygen-and water-free conditions (Figure 3b). Thus, despite a partial oxidation and modification of their magnetic properties, ICNPs maintain high SARs values after 4 months in water and under air, making them promising candidates for biomedical applications in which magnetic heating is required. ICNPs functionalized by the cationic ligand DOP-TEG-C6 display an acute toxicity even for short incubation times. This toxicity can be attributed to the C6 hydrophobic chain carried by the quaternary ammonium group. In sharp contrast, the ligands DOP-TEG-COOH and DOP-TEG-Zwitter provide ICNPs with a very good biocompatibility. No significant cytotoxicity was observed for incubation times up to 4 h at any of the concentrations in NPs of the incubation media. For 24 h incubation, the biocompatibility of ICNPs is very good up to 0.025 mg•mL -1 , the toxicity increasing slightly for higher concentrations. These values are comparable to what can be observed with typical biocompatible iron oxide nanoparticles. [START_REF] Hachani | Functionalisation, and Biocompatibility Studies of Superparamagnetic Iron Oxide Nanoparticles As Potential MRI Contrast Agents[END_REF] Concerning iron carbide nanoparticles, studies are scarce, and until now limited to Fe5C2 nanoparticles. For example, Xie et al observed a 20% cellular death after 24 h incubation with 0.025 mg•mL -1 of polydisperse polymer coated Fe5C2 nanoparticles. [START_REF] Tang | Fe5C2 Nanoparticles with High MRI Contrast Enhancement for Tumor Imaging[END_REF] With similar nanoparticles coated with BSA, Hou et al evidenced the death of 10% of the cells after 24 h incubation with a NPs concentration of 0.025 mg•mL -1 . [START_REF] Yu | Multistimuli-Regulated Photochemothermal Cancer Therapy Remotely Controlled via Fe5C2 Nanoparticles[END_REF] The biocompatibility of our system is significantly better, since no cytotoxicity was evidenced when cells were incubated for 24 h with 0.025 mg•mL -1 Fe2.2C ICNPs functionalized with DOP-TEG-COOH and Zwitter. These results clearly evidence the critical role played by the ligands in the final biocompatibility of the couple [ICNPs-Ligand]. The three systems were further characterized through a cellular uptake study.

HeLa cells were incubated for 4 h in media containing different concentration of functionalized

ICNPs. After incubation, the cells were washed with PBS, and subjected to ICP analysis to determine their iron content (Table 2). 2.0 93 234 [13] 60 [12] For the ICNP-DOP-TEG-C6 system, the concentration of iron inside the cells increases with the increase of the concentration, evidencing the ability of this system to easily enter in the cells. The high cellular-uptake of ICNPs-DOP-TEG-C6 is not surprising, since NPs with positively charged surface interact strongly with the negatively charged surface of cells. [START_REF] Kralj | Effect of Surface Charge on the Cellular Uptake of Fluorescent Magnetic Nanoparticles[END_REF][START_REF] Verma | Effect of Surface Properties on Nanoparticle-Cell Interactions[END_REF] The ICNPs-DOP-TEG-Zwitter system displays a similar behavior, though a drop in the cellular uptake can be observed for the most concentrated incubation medium (2 mg•mL -1 ). This decrease is attributed to a partial aggregation and precipitation of ICNPs in the incubation medium at this concentration. The cellular-uptake of ICNPs functionalized by DOP-TEG-C6

Table 2. Cellular-uptake of different [ICNPs-ligand] systems as a function of the concentration of the incubation media (incubation time = 4 h).

and Zwitter is in the same range than what is usually observed with conventional iron oxide nanoparticles. [START_REF] Riegler | Superparamagnetic Iron Oxide Nanoparticle Targeting of MSCs in Vascular Injury[END_REF][START_REF] Blanco-Andujar | Real-Time Tracking of Delayed-Onset Cellular Apoptosis Induced by Intracellular Magnetic Hyperthermia[END_REF] In strong contrast, no significant cellular uptake was observed for ICNPs functionalized by DOP-TEG-COOH. This result is expected since at the pH at which this study was performed (7.2), the ligands are under the carboxylate form (anionic). Both ICNPs surface and cells surface are thus negatively charged, and the repulsive electrostatic interactions between them prevent the penetration of the NPs in the cells. These observations prove once again that the behavior of nanoparticles in biological conditions is highly impacted by the ligands that are present at their surface.

ICNPs functionalized with the DOP-TEG-Zwitter ligand appear as promising candidates for medical magnetic hyperthermia since they combine high SAR in water, low cytotoxicity and good cellular uptake. However, the environment in cells and tumors is known to have a higher viscosity than pure water, which can affect the heating properties of magnetic nanoparticles. To investigate the influence of the medium's viscosity on the heating power of ICNPs, their SAR was measured in aqueous solutions containing different concentrations of glycerol (Figure 5). The results show that the SAR of ICNPs is slowly decreasing when the proportion of glycerol is increasing, evidencing a dependence of the SAR on the solvent's viscosity. Nevertheless, the SAR values remain high (exceeding 515 W/g up to 20 wt% glycerol, and 320 W/g at 50 wt% glycerol (47 mT, 100 kHz)), making these NPs of interest for practical applications.

To finish, ICNPs functionalized by the DOP-TEG-Zwitter ligand were further characterized by high resolution TEM and STEM-HAADF-EDX at t0, and after 40 days in water (Figure 6). The results show in both cases a core-shell structure with an oxygen rich shell and an iron rich core, suggesting the formation of a thin shell of iron oxide (ca. 2 nm) around the iron carbide core. The presence of this oxide shell already at t0 shows that the surface oxidation of ICNPs is very fast when they are exposed to air and water during the ligand exchange reaction.

Interestingly, no significant difference can be observed between the ICNPs at t0 and after 40 days in water, indicating that ICNPs can be conserved a long time in aqueous solution without suffering from advanced oxidation and losing their advantageous magnetic properties. In addition, their ability to undergo surface oxidation suggest that they could be slowly degraded by liver cells, which is an important factor when considering practical in-vivo applications. [START_REF] Gu | In Vivo Clearance and Toxicity of Monodisperse Iron Oxide Nanocrystals[END_REF][START_REF] Feng | Distribution, Clearance, and Toxicity of Iron Oxide Nanoparticles with Different Sizes and Coatings[END_REF]

CONCLUSIONS

In summary, we show that engineered dopamine-based ligands are suitable for the transfer of magnetic iron carbide nanoparticles into water. The functionalized ICNPs are highly soluble and stable several months in all types of aqueous media investigated (water, PB, PBS, Serum).

Despite partial oxidation, ICNPs maintain good magnetic properties, and most importantly high SAR in water and water-glycerol mixtures. The nature of the ligands is a key parameter governing the cytotoxicity and cellular-uptake of the functionalized nanoparticles. Through proper choice of ligand structure we were in particular able to obtain [ICNPs-ligand] systems displaying very low cytotoxicity, and entering easily or not in the cells.

The biocompatibility of these ICNPs is comparable to the best of the IONs, with their superior magnetic properties making them promising candidates for biomedical applications, and in particular medical magnetic hyperthermia.

METHODS/EXPERIMENTAL

General

All syntheses of iron-based nanoparticles were performed under argon either by using Schlenk techniques or in a glove box. Mesitylene (99%), toluene (99%), and tetrahydrofuran (99%) were purchased from VWR Prolabo, then purified on alumina desiccant and degassed through three freeze-pump-thaw cycles. The commercial products, hexadecylamine (HDA, 99%) and palmitic acid (PA, 99%) were purchased from Sigma-Aldrich. These compounds were used without any additional purification. Dopamine-based ligands were synthesized following methods reported in the literature. Ligand exchange reactions were performed using methanol and dichloromethane degassed by nitrogen bubbling.

Characterization

The size and the morphology of as-synthesized nanoparticles were studied by transmission electronic microscopy TEM. TEM grids were prepared by drop casting a colloidal solution containing as prepared NPs on a copper grid covered with amorphous carbon. Conventional bright-field images were performed using JEOL microscopes (Model 1400 and 2010) working

Hyperthermia measurements

For a typical hyperthermia experiment, an air-tight tube containing about 10 mg of powder of iron carbide or iron/iron carbide nanocrystals dispersed in 0.5 ml of mesitylene was filled under inert atmosphere. The tube was then placed in a calorimeter containing 2.5 ml of deionized water, the temperature of which was monitored during the experiment. The calorimeter was exposed to an alternative magnetic field for a time varying between 10 and 40 s so that the temperature rise never exceeded 20°C. The temperature rise at the end of the magnetic field application was always measured after shaking the calorimeter to ensure the temperature homogeneity, which was measured by two probes (at the top and the bottom of the calorimeter).

The temperature rise was determined after this process from the mean slope of the ΔT/Δt function. Then the raw SAR values were calculated using the expression:

𝑆𝐴𝑅 = ∑ C pi m i i m Fe × ΔT Δt
where Cpi and mi are the specific heat capacity and the mass for each component respectively (Cp = 449 J kg -1 K -1 for Fe NPs, Cp = 1750 J kg -1 K -1 for mesitylene, Cp = 4186 J kg -1 K -1 for water, Cp = 2430 J kg -1 K -1 for glycerol and Cp = 720 J kg -1 K -1 for glass), and mFe is the mass of the pure iron carbide nanocristals.

Synthesis of Fe(0) nanoparticles

In a typical synthesis, the Fe(0) NPs were prepared as follows: in the glove box, 0.65 mmol of PA (333.2 mg) and 0.5 mmol of HDA (241.

Synthesis of iron carbide nanoparticles (ICNPs)

Iron carbide NPs were obtained through the carbidization of preformed Fe(0) NPs. In a typical experiment, 12.5 nm Fe(0) NPs (50 mg, 0.45 mmol of iron) were dispersed in mesitylene (9 mL), and the mixture was pressurized with CO / H2 (2 bar / 2 bar) at 150°C for 140 h. At the end of the reaction, the NPs were recovered by decantation assisted by a magnet and were washed 3 times with toluene (3x5 mL). The NPs were then dried under vacuum. According to the elemental analysis, the black powder obtained contains ~70 wt% of iron. The NPs were further characterized by XRD, VSM and TEM.

Ligand exchange reaction

In a typical ligand exchange reaction, ICNPs were dispersed in a mixture of dichloromethane and methanol (5/1). The water-soluble ligand (5 mass equivalents) was added to the dispersion, followed by triethylamine (250 µL) addition. The mixture was stirred overnight at room temperature. The work up of the reaction consisted of 3 washing cycles with hexane, followed by redispersion in water and purification by dialysis during 48 h. The presence of the desired ligands at the surface of the ICNPs was confirmed by LDI/MS experiments.

Cell viability assay

Cell viability (HeLa cells) was measured by quantifying the reduction of a dye indicator alamar blue. Alamar blue is a dye that takes advantage of mitochondrial reductases to change from oxidized indigo blue state to reduced pink state. This dye has been successfully used in various cell cultures to measure cell viability. Briefly, cells in 4 well plates are loaded with 300 mL culture medium containing 10% alamar blue (Biosource International) and incubated in 37°C, 5% CO2 for 2 h. 100 mL medium from each well was then transferred to a 96 well plate and subjected to measurement. The reduction of alamar blue was measured and calculated by a SpectroMax M5 micro-plate reader (Molecular Device) at 570 nm and 600 nm wavelengths.

₴Max Planck Institute for Chemical Energy Conversion, Stiftstraße 34-36, 45470 Mülheim an der Ruhr, Germany.

Figure 1 :

 1 Figure 1: Dopamine-based ligand structural design.

Table 1 .

 1 Summary of the ligands used to transfer ICNPs in aqueous media. a PB = Phosphate Buffer 5 mM; b PBS = Phosphate-Buffered Saline 150 mM; c Serum containing 10% Bovine Serum Albumin. (I) DOP-TEG-OH; (II) DOP-TEG-COOH; (III) DOP-TEG-TTMA; (IV) DOP-TEG-C6; (V) DOP-TEG-Zwitter.

Figure 2 :

 2 Figure 2: (a,b) ICNPs in mesitylene, (c,d) functionalized ICNPs (DOP-TEG-C6 ligand) in

Figure 3 :

 3 Figure 3: Magnetic properties of ICNPs before (as synthesized) and after 4 months exposure

  Having a library of ligands displaying different functions and charges allows tuning of the cytotoxicity and cellular-uptake of the [ICNPs-Ligand] couple. The cytotoxicity against HeLa cells of ICNPs functionalized with the three different ligands was investigated for several incubation times and ICNPs concentrations ranging from 0.01 to 2 mg•mL -1 (Figure4).

Figure 4 :

 4 Figure 4: Cytotoxicity of different [ICNPs-ligand] systems as a function of the incubation time

Figure 5 .

 5 Figure 5. SAR measurement at 100 kHz for ICNPs (DOP-TEG-Zwitter ligand) dispersed in

Figure 6 .

 6 Figure 6. High resolution electron microscopy of ICNPs in water: at t0 after functionalization,

  at 120 kV and 200 kV respectively. HR-TEM, STEM-HAADF and EDX measurements were carried out using a JEOL JEM-ARM200F Cold FEG working at 200 kV. XRD measurements were performed on a PANalytical Empyrean diffractometer using Co-Kα radiation (λ=0.1789 nm) at 45 kV and 40 mA. Magnetic measurements were performed on a Vibrating Sample Magnetometer (VSM, Quantum Device PPMS Evercool II). XRD and VSM studies were carried out on compact powder samples as well as on diluted assemblies that were prepared and sealed under argon atmosphere. DLS analyses were performed on a Malvern Zetasizer Nano ZS equipment ([0,3 nm -10 µm]). Laser He-Ne 633 nm, 4 mW maxi). Synthesized organic ligands were characterized by[START_REF] Fernandes Cardoso | Advances in Magnetic Nanoparticles for Biomedical Applications[END_REF] H NMR (Bruker 400 MHz). LDI-MS were performed on a Bruker Autoflex III MALDI-TOF mass spectrometer equiped with a Smartbeam 2 laser and a MTP 384 steel target. Analysis conditions : ion source 1 = 19.00 kV, ion source 2 = 16.60 kV, lens voltage = 8.44 kV, reflector voltage = 20.00 kV, reflector voltage 2 = 9.69 kV, pulsed ion extraction time = 10 ns, suppression = 100 Da. TGA were carried out with a Mettler Toledo "TGA/DSC 1 STAR System".

  5 mg) were added to a green solution of 0.5 mmol {Fe[N(SiMe3)2]2}2 (376.5 mg) in 20 ml of distilled and degassed mesitylene in a Fischer Porter bottle. The bottle was then pressurized with dihydrogen (2 bar) and placed in an oil bath at 150°C for 48 h under magnetic stirring (400 rpm). After 48 h, the reaction was stopped and the NPs were recovered by decantation assisted by a magnet, and washed 3 times (3*10 mL) with toluene and 3 times (3*10 mL) with THF. The NPs were then dried under vacuum. According to elemental analysis, the black powder obtained contains ~50 wt% of iron. The NPs were further characterized by XRD, VSM and TEM.

ACKNOWLEDGMENTS

The authors thank ERC Advanced Grant (NANOSONWINGS 2009-246763) and ERC Advanced Grant (MONACAT 2015-694159) for financial support. VR acknowledges the National institutes of Health (US) (EB022641). SUPPORTING INFORMATION. Detailed methods and experimental section, additional characterization of the nanoparticles and ligands.

AUTHOR INFORMATION

Present Addresses

Highly magnetic Fe2.2C iron carbide nanoparticles displaying excellent specific absorption rates were successfully stabilized in aqueous media using specific dopamine-based ligands. Even after several months of exposure to water and air, iron carbide nanoparticles retained high heating power through magnetic induction. The nature of the ligands was proven to be a critical parameter to tune the biocompatibility and cellular-uptake of the considered nanoparticles.


X Liu

B Revol

F Costa

PARASITIC ELEMENTS MODELING AND EXPERIMENTAL IDENTIFICATION IN A GAN HEMT BASED POWER MODULE

In this paper, the parasitic elements of a GaN based power module are investigated. Because the very high-speed switching of GaN high electron mobility transistors leads to overvoltage, conducted and radiated emissions, it is crucial to manage accurately the routing of the power module in order to minimize parasitic capacitances and inductances. The paper focuses on the prediction of these parasitic effects and their comparisons to measurements. In this approach, the difficulty lies in the determination of very small values of the parasitic elements which requires specific operating for calibration and measurements.

I. INTRODUCTION

The emergence of wide bandgap GaN high electron mobility transistors (HEMT) is promising for power electronics performances if their packaging is matched to their very high switching speed [1], [START_REF] Reusch | Understanding the Effect of PCB Layout on Circuit Performance in a High-Frequency Gallium-Nitride-Based Point of Load Converter[END_REF]. Actually, only discrete components can be found out of-the-shelve with low-inductance packaging (PQFN or similar), but no power module are already available. Our purpose is to study such a new type of module exhibiting low parasitic elements and including capacitors for HF differential and common mode decoupling, in order to get the best performances as possible and to limit conducted high frequency emissions [START_REF] Wolfgang | A race towards the lowest inductive power module[END_REF]. For this objective, we have realized a prototype of power module using discrete HEMT's as dies (GaN System GS66508P) in order to develop our simulation and characterization methodologies. The first part of the paper will describe the bare module and how to calculate and measure the parasitic capacitances and inductances. In particular, it will be shown how to extract the module parasitic capacitances even if the dies are soldered in. Comparisons will be presented that validate the methodology. The second part will present the electrical model of the module including decoupling ceramic capacitors and how it is used for EMI simulation of the CM current issued of the module.

II. BARE MODULE MODELING

A bare module of a GaN-based switching cell was designed in the lab with neither in-module multilayer ceramic capacitors (MLCC) nor semiconductor devices.

Its layout is presented in Fig. 1, it is modeled as the circuit in Fig. 2. C M+ , C M-and C MA are the main parasitic capacitances and L M+ , L M-are the equivalent parasitic inductances of the copper tracks of the DC bus. The parasitic elements can be either calculated by a numerical approach (3D PEEC method) such as ANSYS Q3D (c.f. II.1-A) or be measured directly by an impedance analyzer with proper protocol. It should be noted that both methods lead to some error. However, results-matching between the two methods makes a tendency to a good estimation. Last but not the least, their values will be finally evaluated by comparison between circuit simulation and measurement in a real test bench. The capacitances between different copper pads on top layer are too tiny and therefore they will be ignored. It should be noted that C SS (capacitance between the kelvin source pad and GND) should be added with C M-and C MA respectively since they are connected inside the power device's package.

Tab. 1: Q3D extraction of parasitic capacitances in our bare module (unit: pF)

C M+ =38.2 C M-=41.7 C MA =34.9 C G =2.5 C SS =2.0 b -Impedance measurement
The parasitic capacitances can be directly measured one by one in the bare module without any semiconductor power chip. However, the bare module (only layout and interconnections, no semiconductors) is not often achievable. The measurement approach below provides the capability to measure parasitic capacitances directly in a completed module including semiconductors. In a first place the two GaN devices should be short-circuited between the gate and the HEMT kelvin source so that each device is equivalent to a capacitance between drain and source with a value of C OSS (C DS +C GD ) at low V DS voltage (~ 500 pF for GS66508P). Therefore there are 5 unknown (C M+ , C M-and C MA ) or roughly known (C OSS of high-side and low-side devices) capacitances between the DC+, DC-, AC and GND potentials (Fig. 3).

Fig. 3: Capacitances between different terminals.

Since the error of several pF between the C OSS of highside device and the low-side one (mostly due to the manufacturing process) can make a difference compared to the values of C M+ , C M-and C MA (according to the Q3D calculations), it is advised for precision purpose to treat the two C OSS as two unknown coefficients -C OSS_H and C OSS_L . Finally 5 capacitance measurements with proper short-or open-circuit operations should be done in order to determine the 5 unknown capacitances between terminals, as discussed hereafter.

In order to measure the capacitances of the module mentioned above without any soldering, a specific board has been realized. It consists of 2×1 spring loaded connector arrays used to take contacts on different terminal areas on the top layer of the module. Each pair of those connectors allows to measure or to short-circuit the capacitance between two terminals. All contact points are shown in Fig. 4a. Each pair of the spring connectors is in parallel with a SMB adaptor for impedance measurement and a 2×1 pin array for short-circuits convenience by a jumper. All the connectors are configured in a PCB which is fixed above the power module via four screws (Fig. 4b). It should be noted that COM contact points (Fig. 4a) are used for the measurement of the stray capacitance in the spring connector together with the SMB connector, which should be compensated (subtracted) from the measurement results.

II.2. Preparation & Measurement method choosing

The measurement procedure is presented hereunder. The contact points CP5 (gate and kelvin source for low-side GaN device) and CP6 (gate and kelvin source for highside GaN device) are short-circuited before all SMB Spring connector Pin array measurement; in addition, since the "to be measured" capacitances' value is hundreds of pF, corresponding to very high impedance at the frequency range up to 30MHz, the "GP series" impedance measurement configuration (Gain-phase with series-through) will be used with our network analyzer (Keysight 5601).

a -Calibration & Compensation:

After O/S/L calibration to the end of the cable, at least an "OPEN" compensation should be done at the COM contact point, in order to eliminate the test fixture stray capacitance between two spring connectors and the SMB connector. The value is measured as 2.0 pF. The final module's model together with the DC link capacitor is presented in Fig. 7. There is an extremely large power loop inductance (over 60 nH) in which the bulk capacitor and the power module terminals (including the capacitor supporting PCB board and the interconnection screws) contribute respectively around 35~40 nH and 22 nH.

III. IN-MODULE CAPACITORS MODELING

As previously exposed, it may be efficient to include a first stage of HF decoupling directly inside the module in order to limit overvoltage on the DC bus and to partially filter CM currents. This CM and DM filtering is shown in fig. 8a

III.2. Final model of the module with in-module ceramic capacitors

At last, the full circuit model of the module including C X and C Y capacitors with DC link bulk capacitor is shown in Fig. 11, with parasitic elements in blue color and discrete components in black color. This full model can be used now to simulate any kind of electrical quantities regarding the working of the module and/or EMI aspects [START_REF] Liu | Methodology for EMC Analysis in a GaN Based Power Module[END_REF]. A specific behavioral model of the GaN components was previously developed in [START_REF] Liu | Methodology for EMC Analysis in a GaN Based Power Module[END_REF].

For example, Fig. 12a presents a comparison of the simulated and measured spectra of the CM current issued of the module without any decoupling capacitor. These curves show good agreements between the predictions and the measurements in a very large bandwidth up to 100 MHz which is the limit fixed for circuit modeling of this study. Over this frequency, modeling would require to take into account other effects not taken into account here over like propagation and radiating effects.

IV. CONCLUSION

In this paper, an experimental method to extract the parasitic elements of a GaN based power module has been presented. This specific method is achieved while power dies are already soldered in the module, requiring an adequate calibration setup specifically dedicated to the module under study. The values of the measured elements have been compared to the ones extracted by a numerical tool; good matching's have been observed. Moreover, the electrical model obtained with these techniques enable very accurate simulations of any electrical or EMI quantities issued of the module over a large frequency range. At last, the method presented hereunder can be extended to other circuits different from power modules but including semiconductors.

Fig. 1 :Fig. 2 :II. 1 .

 121 Fig. 1: Layout of our power module (4 x 4 cm), capacitances to be calculated in Q3D.

Fig. 4 :

 4 Contact points for capacitances' measurement (a) and measurement board (b).

Measurement 1 :II. 3 .Fig. 5 :II. 4 .

 1354 Fig. 6: Comparison between the fitted model and impedance measurement between TP+ and TP-.

Fig. 7 :

 7 Fig. 7: Values of main parasitic elements in bare module with the DC link capacitor.

Fig. 8 :

 8 Fig. 8: Power loop for different C X capacitors (a) and the equivalent circuit (b).

Fig. 9 :Fig. 10 :

 910 Fig. 9: Illustration of the inductive loop of C Y+ /C Y- capacitors.

Fig. 11 :

 11 Fig. 11: Full model of SATIE module with in-module capacitors.

  Fig. 12b shows the impacts on the first spectrum of the C Y decoupling capacitors for two different values. a) b) Fig. 12: Spectra of the module CM current, a) module without any local decoupling, blue: measurement, red: simulation. b) left : C y =0.7 nF, right C y = 7 nF.

ACKNOWLEDGMENTS

This work has been supported by the PIA MEGaN.


Olivier Bui

Xavier Leoncini

Growth of a tree with allocations rules: Part 1 Kinematics

Keywords: 05, 45, -a Nonlinear dynamics and chaos -05, 65, +b Self-organized systems

A non-local model describing the growth of a tree-like transportation network with given allocation rules is proposed. In this model we focus on tree like networks, and the network transports the very resource it needs to build itself. Some general results are given on the viability tree-like networks that produce an amount of resource based on its amount of leaves while having a maintenance cost for each node. Some analytical studies and numerical surveys of the model in "simple" situations are made. The different outcomes are discussed and possible extensions of the model are then discussed.

Introduction

Systems of transportation frequently appear in physics, engineering, biology: transportation of water, electricity and gas in cities, street network, river basins [1] or vascular systems [2,3] etc. A typical way to formalize systems of transportation is to transform it in a problem of flow in a network between sources and sinks [START_REF] Ford | Flows in Networks[END_REF][START_REF] Carmi | [END_REF][6] with discussions about finding the best path from sources to sinks. Such approach can lead to application in economic context [7][8][START_REF] Hougaard | Sharing in Networks[END_REF]. More physical approaches involving transport of fluids or electricity have been considered too and, in these cases, use of the formalism of electrical circuit (Kirchoff's law, current, potential etc.) is not unusual [START_REF] Bohn | [END_REF]11], the more geometrical aspects would also be involved [12,13]. This approach is often focused one minimizing some dissipative energy which gives it some application in engineering context but also in the description of nature as natural networks may arise from minimization of some dissipation [14,15].

We can notice that the focus is typically on optimization. The substance being transported is ether not a resource to be used or consumed or its use is of little importance to the problem the papers are studying. A possible problem with such approach is that real life networks may not just form from minimization/maximization problems but from some interaction between to flow and the network: vascular system of plants is a network transporting sugars, but these sugars are also used to build the network itself. Perennial plants (including trees) do not rely on a centralized organs such as a heart or a brain to grow or distribute the resource throughout the plants and yet how they branch and their final shape is not totally planned at birth either but is also strongly caused by adaptation to external stimuli [16]. Given this, we may believe trees a Present address: Insert the address here if needed are an example of self-organization arising from very simple local rules and interactions [17]. The vascular system of trees is composed of a part called xylem and another one called the phloem. The xylem is mainly responsible for water transportation. This flow of water is unidirectional: it goes from the roots to the leaves then most of it gets evaporated at the level of leaves. These leaves then play a major role for the growth and sustenance of the organism: they are the ones absorbing carbon in the air and create sugars through photosynthesis. These sugars are essential as they sustain the living cells of the tree as well as provide its building materials. The part of the vascular system responsible for transporting this vital resource is the phloem which is the part we are interested with since it is the one that potentially showcases this idea of interaction between the network and the resource being transported. The biophysics of the phloem has been extensively studied [START_REF] Munch | Die stoffbewegungen in der pflanze[END_REF][START_REF] Schepper | [END_REF][20].

However, instead of studying the biomechanics of sugar transport like the aforementioned papers, our approach is more abstract and in more line with the work about flow, sources and sinks mentioned earlier but we also propose to consider as well how the flow can be used to built the network itself while in turn the flow depends on the network and make a simple dynamical system out of it. In order to start on relatively simple grounds the toy-model we built and present in this paper will only use tree network (i.e. loopless). This has both the advantage of keeping the network topology quite tractable, while potentially shed some light on actual possible growth mechanisms occurring in biological trees. The dynamical and growing aspects we added also makes it similar to cellular automata which was a concept introduced by John Von Neumann in [START_REF] Neumann | Theory of self-reproducing automata[END_REF], though a more famous example of cellular automaton is the Game of Life invented by John Conway [START_REF] Gardner | [END_REF]. The tree has 5 extremities and is of height 3. And if the tree follows Leonardo's rule then the volume of T will be 5 units of volume, A will be 2 units and B will be 3 units.

In short, we are studying growing systems for which the growth and self-organization is driven by short range interaction (exchange of resource). The paper is organized as follows, in Sec. 2 we discuss the ingredients of model inspired from trees, trying to keep a minimal set of ingredients or variables. Then in the following section 3, we propose a way and implement the self-consistent dynamics of the growth of the tree. In Sec. 4 and Sec. 5 we perform some analytical calculations and predictions of the model and then perform some more analysis with numerical simulations in Sec. 6 before concluding in Sec. 7.

A bio-inspired model

In this section we construct a "growing" system centered around resource distribution and allocation inspired by biological trees. Considering how complex trees are, we opted for what appeared to be a simpler description of their growth and resource allocation. The main constituent of a tree is carbon, it is absorbed by their leaves and their source of energy is sunlight which also involve leaves. Both the carbon and energy is distributed and allocated in the form of sugars through the phloem. So for a simplistic description of trees we may shave off the xylem, water transport and the roots, which leaves us with a model featuring leaves as our sources, sugars as our resource, branches representing the nodes of our growing transportation network, and the volume of carbon each branch has fixated. Each branch will need to regularly consume the resource (sugar) to increase in volume (of carbon) and to keep its cells functioning. Let us introduce some terminology for our model:

-Trees are modeled as acyclic (i.e. loopless) connected networks like in figure 1, the nodes are called branches, one of them is designated as the trunk which is the only node existing at the start of the simulation.

-Each branch has a height which is an integer indicating its distance to the trunk (the trunk is of height 0, the branches directly linked to it are of height 1 etc). A branch with greater height than another is said to be higher.

-If two branches are linked by an edge then the highest one is said to be the child of the second. The second is the parent of the first. The terms ancestors/descendants will be used for parent/children, grandparent/grandchildren, etc. And we say a branch A descends from another one B when the branch A is a descendant of B. -Finally, we call extremities branches without any children. And the height of the tree is defined as the height of the highest extremity. Now, we propose a few rules related to resource distribution and allocation:

-The extremities will be the sources (the only branches possessing leaves) and all the branches, including the extremities, need to consume sugars as maintenance cost.

-Each branch can also use up its resource to increase its volume, the idea being large volume equates with more resilient branch, but in exchange the maintenance cost is higher for a larger branch. -Each branch can transfer resource to any of its children or to its parent and it can also store it in its "personal" reserve.

And the rules driving changes in the network are the following:

-If a branch dies (could not pay its maintenance cost) then all its descendants will die (the branch is cut off the tree). Since the sources are at the tips of the tree, it creates an interdependence between the ones producing resource and the rest of the tree. -About growth: extremities can use up resource to create children and in exchange cease to be sources while the newborns will be the new extremities. It is through this process that the tree can increase its number of sources since each time a branch creates more than one child, the number of sources increases.

Time is discrete: the branches act in successive rounds. We will provide later an example of how this "game" may progress.

The last characteristic we need to detail is the volume of branches. Each branch possess their own quantity of volume (i.e. they have a certain size). The bigger a branch is the more costly its maintenance is, and to become bigger the branch also needs to spend resources. In our model, we decide that the minimum volume (in some arbitrary unit of volume) a branch needs to have is equal to the number of its descendants that are extremities: if a branch has one child that is an extremity, two grandchildren that are extremities and one great grandchild that is an extremity, then the branch at least needs a volume equal to 4 (in unit of volume). This rule will be called Leonardo's rule because it is inspired from the real Leonardo's rule which has been studied in plants [START_REF] Vinci | The Notebooks of Leonardo da Vinci[END_REF][START_REF] Shinozaki | [END_REF][25][26].

The exact formulation of the real Leonardo's rule is a rule of conservation of cross section areas: when looking Fig. 2. The dynamics follows this sequence of events: the sources get an amount of resource. Then they can choose to store some amount for themselves and give the rest to their parent. Then the parents do the same to the grandparents etc, until the flux reaches the trunk. After this, the trunk pay its maintenance cost and use some of the resource to grow in volume, then divides the remaining to its children. The process is repeated until we reach the extremities with no child, these last one can use up the remaining resource to birth new branches. During the phase where the flux goes upward called "flux up", branches that cannot pay their maintenance or grow to the correct volume may die. at a branch that branches out into several others ones, the cross section of the branch before the branching node is equal to the sum of the cross sections of the branches after the node. Our Leonardo's rule would be synonymous with the real one if consider the branches in our model as cylinder with all the same length and the extremities of our tree all have a volume of 1 unit.

Example of implementation

We have established the basic rules describing our model. However such description does not specify the model we implement. We still have to describe how the "game" progresses in time or, in other words, the dynamics of the system.

Example of a possible dynamics: we start with an amount of resource created at each source then during each "turn" every branches can store, consume for volume growth or transfer any amount of resource to their neighbors. And then after x number of "turns" each branch has to pay maintenance from their stored resource or die and we repeat process. Creating new branches may involve consuming resource for y numbers of unit of time etc. Such dynamics would be fairly difficult to study, so we opted for a simpler dynamics. Unlike the example described above, in our chosen version, the branches do not act simultaneously but in succession. We will call this succession of events and actions a generation. When a generation ends, a new one begins, so generations will be used as our unit of time. Each generation is composed of two phases and each branch acts once during each phase. The first phase describes how the amount of resource generated at the extremities flows "downward": from the extremal branches to eventually the trunk. And the second phase describes how the resource reaching the trunk bounces back "upward": from the trunk to the extremities.

Let us specify this first phase which we call "flux down" (Fig. 2): the first branches to act are the extremities, then the next to act are the branches for which all the children have already acted, and we repeat until every branches have acted. Each action is the following sequence of events: -If the branch is an extremity then it receives an amount p 0 of resource. Otherwise, it receives an amount from its children and will remember which amount was received from which child. -It also receives from each child information about their energetic needs which is their maintenance cost as well as the resource they need to create the additional volume of wood necessary to support the structure above them. This information about energetic needs determines how the resource will be distributed during the second phase when the parent will have to choose the amount to distribute to each child. -The branch can keep some of the resource in its reserve then transfer the rest to its parent. Then, it is this parent that will perform this same sequence of events. In the second phase which we call "flux up", the branches also act in succession but the order is reversed (Fig. 2) and the actions are now:

-If the branch is the trunk then it will start with all the resource that was given to it during flux down, then the flux will bounce back. If it is not the trunk it will have the resource it kept in its reserve but has also received an amount from its parent. -With this amount of resource, it must pay its maintenance cost which depends of its volume and dies if it can not. -Then, it will try to use up resource to grow to its "optimal" volume which is equal to the number of extremities that are its descendants (the volume of a newborn branch being set to one unit of volume, it means if the volume of each branch is "optimal", in the sense we have defined above, then the tree will perfectly respect Leonardo's rule). If the branch volume is not larger or equal to its optimal volume, it is cut off. -The branch shares all its remaining resource to its children. The share given to each one depends on what happened during the flux down phase: the energetic needs each one told as well as the amount of resource each gave. -If the branch has no child (it is a source) then it uses the resource to create children. There is a maximum number of children it can create. In this phase, the proportion a parent shares/distributes its resource to one child or another is decided by a formula of our choosing. However the parent is supposed to remember the amount of resource each child gave during flux down as well as their energetic need. So a suitable formula would be one that use these two values, doing so we obtain the tree represented in Fig. 3.

In the end the dynamics of the tree only depends on a few parameters: Using the simple dynamics we described, we can obtain the tree that was represented drawn in three dimensions in Fig. 3. All the branches were identical: during flux down, they all use the incoming resource using the same algorithm (they grow in volume in order to get the volume dictated by Leonardo's rule) during flux up, they all use the same algorithm to decide how they share their resource among the children (the algorithm take as entries the energetic need of a child and its previous contribution during flux down) but with a very small proportion of the resource being sometimes distributed randomly. As for the branches with no child, they would use the remaining resource to create as many children as possible.

On the other hand, if, during flux up, we do not implement that slight random distribution then the branches will always split equally their resource among their children because every branch have the same strategy so, during flux up redistribution of resource to the children, the parent will have no reason to discriminate among the children. This results in "symmetric" trees. Either way, the growth and trees we obtain depend on the parameters (p 0 , C, m 0 , C v , α and N max ). Therefore, we may want to have a theoretical view on the system. The simplest way to look at the system is by having a static approach: for instance, for a tree with a given topology we can calculate the total resource produced and its total maintenance cost assuming the volume of the branches follows Leonardo's rule.

A static description of the problem for the case α = 1

If we solely look at the total production and maintenance cost of a tree for a given topology, it will allow us to make statements that are independent of the dynamics we chose, as such, this is the approach we will start with. For the calculations in this Section, the only hypothesis that will be used are: every extremities produce the fixed amount p 0 per unit of time, each branch pays the cost m 0 × (Vol branch) α per unit of time and the volume of the branches are such that they respect Leonardo's rule, i.e no branches have been cut above it.

We may want to interpret α. If we were look at our tree in the context of botany: the branches in our model are supposed to be cylinders of same length and the only morphological difference between the branches is their radius, therefore the volume scales like its cross-section surface. α = 1 means the maintenance cost is m 0 × (Vol branch): this proportionality to the volume can be interpreted as the cells being uniformly spread in the branch and consuming resource at the same rate. In reality, for large trees and branches, living cells are thinly located at the exterior of the branch while the interior of the branch is mostly dead wood: this would correspond to α 1/2. Outside the context of botany an analogy to be drawn between our cost maintaining a structure (volume) and the results on dissipative cost of transporting resource found in [27]. (The comparison mainly holds thanks to the conservation of volume i.e. Leonardo's rule.) Their paper found trees (loopless network) as the "optimal" topology when 0 < α < 1 but not in other cases. Considering this fact, we will ignore α > 1. The case limit α = 1, however, can still give us some insight for 0 < α < 1 and makes analytical calculations easier.

Since we assume all extremities (and only them) create an amount of resource p 0 , then the total production, noted P, is equal to p 0 × E where E is the total number of extremities the tree has. Now, we recall that, for each branch, its height is an integer defined as the distance between the branch and the trunk within the graph representing the tree. For each extremity of the tree, we can note its height, then we define E h as the number of extremities that have a height equal to h, as a result E = E h . Thus, by defining H as the height of the highest extremity, we can write:

P = H h=0 E h × p 0 (1)
We call M the total maintenance cost. It has a simple expression when α = 1 (the linear case):

M = H h=0 (h + 1) × E h × m 0 (2)
To illustrate this formula we can look back at Fig. 1: if Leonardo's rule is respected the trunk has a volume of 5 because there are 5 extremities that are its "descendants". Again with Leonardo's rule, if we look at the children of the trunk (the branches named A and B): one has a volume equal to 2 and the other to 3. In other words, summing the volume of the children of A and B also yields 5 which is the number of extremities that either "descends" from A or from B. Then if we look the branches at the height just above, the sum is also 5 which is the number of extremities descending from them plus the number of extremities among them. This last reasoning yields 2 for the branches of height 3. This observation can be generalized into the following rule: for a tree of height H and given any integer j ∈ [0, H], the sum of the volume of all branches that have a distance j from the trunk (height equal to j) is equal to A j = H h=j E h . So the total of volume of the tree is H j=0 A j which can be rewritten as

H h=0 (h + 1) × E h .
A way to interpret this last expression is by saying that adding an extremity of height h increases the total volume by (h + 1) units. Thus, if the maintenance cost of a branch is proportional to its volume (i.e. α = 1) then, by linearity, we get the equation ( 2) by multiplying the total volume of the tree with the proportionality factor m 0 . Hence for α = 1, by combining ( 1) and ( 2) we obtain the global balance B:

B = P -M = H h=0 E h (p 0 -(h + 1)m 0 ) (3)
This expression for B is simple enough we can easily deduce a few results. We recall N max is the maximum number of children a branch may have. For a fixed N max < ∞, and using the expression (3), we obtain the following results:

1. We can not find arbitrarily tall tree such that B ≥ 0.

This in turn implies that no tree can grow infinitely in height as they would reach a height above which B can only ever be strictly negative. 2. For trees of height H ≤ H 0 , with H 0 being an integer we will define later in Eq. ( 4), we know the form of the tree maximizing B ≥ 0: this maximal tree is the tree for which every branches have N max children, except the branches of height H. I.e. this is the tree of height H that has N H max , hence maximizing its number of extremities.

In order to show this we start by defining the integer h max = p0 m0 -1. The brackets [x], indicate the integer part of x, and will be used as such hereafter. This is the smallest integer for h such that p 0 -(h + 1)m 0 is negative for any h > h max . Therefore, in Eq. ( 3), the terms of indexes higher than h max are negative. In other words, any extremities added at a height strictly higher than H max will contribute negatively to our profits B. Consequently, since N max < ∞, it is clear we can not find arbitrarily tall trees such that B as we would reach a point where adding branches will only contribute negatively. Now, we prove the second result: let us consider a tree. We consider a branch such that its height is strictly less than H max and has a number of children n c such that 1 < n c < N max (so we precisely pick a branch that is not an extremity). If we build a tree identical to this first one except that the aforementioned branch is given new children until it has N max children (with the added children being extremities (i.e. have no child) ) then the second tree Fig. 4. α = 1. The tree on the on the right is obtained from one on the left by adding the branch highlighted in red. The equation (3) implies that, by giving more children to non extremal branches of height < Hmax, we will always increase B. In other words, if Hmax > 2 then the tree on the right has a bigger B then the one on the left.

will have a higher B than the first one. An illustration of this statement is in Fig. [START_REF] Ford | Flows in Networks[END_REF]. The second step is to take, this time, an extremity and to give it N children which will become the new extremities, then we can look at how B changes. Since we now have N -1 more extremities, the production increases by (N -1)p 0 but the maintenance cost increases by ((N -1)(h + 1) + N )m 0 . With a bit of algebra we can easily determine the values of h for which the increase in production compensates the increase in maintenance: if the height of the extremity we are giving new children is equal or lower than h 0 = p0 m0 -N N -1 -1, then the contribution to B is positive. It is negative otherwise, and leads to

H 0 = p 0 m 0 - N max N max -1 (4)
By combining the two previous remarks, we conclude that, below a certain height, the tree with the highest B is the tree with N H max extremities and, beyond this height, any adding new children to a extremity lowers B.

Generalization for 0 < α < 1 and infinite trees

To sum up, the two important facts about α = 1 we have established are: no infinitely growing tree is viable and creating as much children as possible yield the most productive trees as long as the height of the tree is under a certain height H 0 . Now, for 0 < α < 1, we may ask whether arbitrarily tall trees that are viable exist. In other words, we want to know whether, for each height H, there exists one tree such that B > 0 and what constraints on p 0 and m 0 are needed. Obviously, calculating B for all possible trees of height H may not be the best strategy. However, inspired from the results of the α = 1 case, we can first limit our search to trees of height H with N H max extremities (trees maximizing its number of branches). To illustrate how B can be calculated, let us take the exampleN max = 3 and H = 3 then we can easily count that the trunk has a volume of 27, its 3 children have a volume of 9, its 9 grandchildren have a volume of 3 and, at last, there are 27 extremities. Therefore, the total maintenance is:

M = m 0 (27 α + 3 × 9 α + 9 × 3 α + 27 × 1 α ) = m 0 3 i=0 3 i × 3 (3-i)α M = 3 3α m 0 3 i=0 3 (1-α)i (5)
Using the same kind of reasoning for a more general N max (that we will note as N from now on) and any H, we can deduce the Eq.␣ (6):

M H = N αH m 0 H i=0 N (1-α)i = N H m 0 N 1-α -N -(1-α)H N 1-α -1 , if N > 1 and α = 1 (6)
From there we have an expression for B H when N > 1 and α = 1 :

B H = N H p 0 -m 0 N 1-α -N -(1-α)H N 1-α -1 (7)
From this expression we can see that if

0 < α < 1, then p 0 -m 0 N 1-α -N -(1-α)H N 1-α -1 ≥ p 0 -m 0 N 1-α N 1-α -1
for every H therefore, for a given N and α, the value of p 0 /m 0 determines whether B H will eventually become negative. So, we have the following two cases:

p critical = N 1-α N 1-α -1 p 0 /m 0 ≥ p critical ⇒ ∀H, B H ≥ 0 p 0 /m 0 < p critical ⇒ ∃H f , ∀H ≥ H f , B H < 0 (8)
In the first case, B H → ∞ and, in the second case, B H becomes strictly decreasing after a certain H, then B H → -∞. So, a sufficient condition for the existence of viable trees (i.e. B > 0) at any height is p 0 /m 0 ≥ p critical .

This condition is also necessary, but we need to go through multiple arguments to reach this conclusion. A first observation to make is that branches have no benefit having only one child since only having one do not increase the number of extremities while increasing the maintenance. So, we will only discuss trees for which every branches, except the extremities, have at least 2 children. With this kind of trees, for any m ∈ N, we can find H m ∈ N such that any trees of height H m or higher have m or more extremities.

We make a second observation: let m ∈ N, we can find H ∈ N such that N H ≤ m < N H+1 . If we construct a tree with m extremities by starting from the tree of height H that has N H extremities, then adding the remaining r 0 = m -N H new extremities in the way described in Fig. [START_REF] Carmi | [END_REF], this tree has a maintenance cost lower or equal than any other trees with m extremities. Let us verify this last statement. First, it is obvious we there can not exist a tree with m extremities that is strictly shorter than the one we just constructed because no branch may have more than N Fig. 5. Here, N,the maximal number of children a branch may have, is equal to 3. If we want to construct a tree with m extremities (in this example m = 18), then one way it can be done is by starting from tree of height H (in our example, H = 3) maximizing its number of children then we give children to consecutive extremities until we reach m, like in the figure above. This way of constructing a tree with m extremities minimizes the total amount of volume the tree needs to have (assuming Leonardo's rule applies).

children. Second, let us consider strictly taller trees with m extremities: they would have more total volume (since we recall the formula for the total volume is h (h+1)E h ) and share this total volume among more branches, therefore, using the inequality ( x) α ≤

x α , we deduce their maintenance cost must be higher. Third, now we simply need to compare trees which height is H + 1: they all look like our tree except that the r 0 remaining extremities may be split differently. Our splitting minimizes maintenance costs because v α + (v + x 1 + x 2 ) α ≤ (v + x 1 ) α + (v + x 2 ) α which means it is better to pack the r 0 branches together in one side.

A third and last observation has to be made. We mentioned that when the condition is not respected (i.e. p 0 /m 0 < p critical ), (B H ) H is strictly decreasing and negative after some H f . We can establish that if H ≥ H f , then the tree with m (such that N H ≤ m < N H+1 ) extremities described earlier (Fig. ( 5)) must have a B inferior to B H (therefore negative too). Indeed, B not being inferior to B H iscontradictory: we name again r 0 = m -N H and assume r 0 = N -1 then conclude that if the resulting B is superior to B H , then B H+1 ≥ B H because in this case adding one (or more) group of N -1 extremities will increase B even more. Consequently, the case r 0 = N -1 must lowers B which we can then use to prove the case r 0 = k(N -1) also lowers it, and after this we deduce it for a more general r 0 .

Finally, using the first observation, we can find H 1 such that every trees of height higher than H 1 have at least N H0 extremities. Then, for each of these tree, using the second observation, we can create a tree, similar to the one in Fig. [START_REF] Carmi | [END_REF] with the same number of extremities but with a lower B. But, this last tree has a height H ≥ H 0 so its B is negative, using the third observation. Thus, when p 0 /m 0 < p critical , arbitrarily tall and viable trees do not exist. Now, let us go back to Eq. ( 7) and contextualize it. Until now, we looked at the static cost of the branches. By doing so we ignored an important aspect of our dynamical tree: branches need to actually spend resource to grow into the volume of Leonardo's rule; each time a branch increases in volume it consumes an amount proportional to this volume increase (Sec. 3: C v being the proportionality coefficient). Let us look at the cost induced by this process when we go from the tree of height H -1 and N H-1 extremities to the tree of height H with N H extremities.

C vol(H-1→H) = C v (number of additional extremities)H = N H C v (1 -N -1 )H (9) If we can compare C vol(H-1→H) to B H , it is C vol(H-1→H)
that dominates for large H (B -C vol tends to -∞ as H grows). However unlike a maintenance cost, this creation cost is only paid once each time the tree grows, so a tree that takes time to accumulate a reserve between each growth spurt may grow arbitrarily tall. But the amount of time waiting to accumulate a reserve will increase as O(H). Simply put, if C v = 0 then our tree will not be able to grow infinitely tall unless it is "intelligent" enough to slow down its growth, and even then we predict a "rate of slow down" that is exponential. It may echo the fact real trees have their growth slowing down with size [28].

6 Simulation of symmetric trees

Encoding the trees in sequences of integer

Up until now, all our theoretical results are independent of the dynamical rules we implement: they all are simple arguments on topology and economic costs. Let us go back to the dynamical system described in section 3, and we consider the case where all the branches has the same strategies during flux up (the phase in which the flow go from the trunk to the extremities): at the end of each generation, the extremities generate as much children as possible with the resource they possess. As explained in Sec. 3: without asymmetry or randomness, the tree will be "symmetric", meaning that all branches located at the same height will have the same number of children and brothers, the same volume etc. The trunk will have u 0 ≤ N number of children. By symmetry, these u 0 children will all have the same number of offspring, noted u 1 ≤ N , etc. So, for a tree of height H, we can construct the sequence u 0 ,u 1 ,...,u H-1 where u i ≤ N . A tree composed of only a trunk is represented by an empty sequence. Conversely such a sequence will define a symmetric tree of height H.

Hence we can see the growth of our simulated symmetric trees as a finite sequence of u i ≤ N that evolves dynamically. While the "form" of a symmetric tree of height H is characterized by a sequence u 0 ,...,u H-1 , the tree itself is also defined through how the volume and the reserve is allocated among the branches. Since every branches located at the same height have the same volume and reserve then the full characterization of a tree of total height H is given by 3 finite sequences: (u i ) 0≤i≤H-1 , (v i ) 0≤i≤H representing the volume v i of the branches located at height i and (r i ) 0≤i≤H giving the amount of reserve r i stored in the branches located at height i. With our symmetrical dynamics the sequence (u i ) 0≤i≤H-1 may shrink (branches of height superior to some h die) or may grow (new integers u H , u H+1 etc added) but the values u i themselves do not change because if a branch of height h + 1 dies, those of the same height also all die by symmetry, therefore u h would simply disappear from the sequence as well as all u i>h .

The initial condition for all our simulation is (u i ) = ∅ (symbol for empty sequence), (v i ) = ( 1) and (r i ) = (0) (i.e. we always start with a tree formed only by a trunk of volume 1 with no initial reserve). Given how, we kill every branch that has a volume less than what Leonardo's rule would predict then, as already mentioned, the only time some v i do not follow Leonardo's rule is when a tree got some branches dead (so that the survivors will have bigger volume than predicted by Leonardo's rule). So if we interest ourselves with infinitely and steadily growing trees the sequence (v i ) i is unneeded since it matches Leonardo's rule perfectly. Though, C v = 0 is a requirement if we want trees growing ever steadily.

Evolution of (u i ) i for infinitely and steadily growing trees

We define a "steadily" growing tree as a tree such that at the end of generation number i, the height of the tree is also i (i.e. the tree never shrinks and never stagnates).

We have stated at the end of Sec. 5 that unless C v = 0 an infinitely and "steadily" growing tree is impossible, so from now on C v is taken as 0 (i.e. growing in diameter/volume do not cost any resource to the branches). If we note B H = u 0 × u 1 × ... × u H-1 and define B 0 = 1 then the B i represent the total number of branches located at the height i. If we want to write an expression for u H+1 or simulate the system, we first need to place ourselves in a specific version of the system described in Sec. 3. For the specific system we study: during "flux down", the children give everything to their parent which means that at the beginning of "flux up" the trunk start with the total production. Then, during "flux up", the parents only pay their maintenance, are eliminated if they can not pay and then grow in volume (which is free since C v = 0) but do not keep anything more in reserve and distribute equally to their children the remaining flux. Finally when the flux reaches the extremities, they, after paying their maintenance, use as much resource as possible to create as many children as they can (the cost of creation for each child being an integer C). With this choice of evolution rule, the system is greatly simplified because, right before the extremities have to create new children and at the end of all the maintenance payment, the resource the extremities have is equal to the total amount produced (by the extremities at the beginning of "flux down") minus the total maintenance cost of the tree plus the leftover resource that could not be transformed into children the turn before, noted R i . So since we know the resource the extremities have right before creating children, we can predict the number of children each extremity will spawn, u H+1 :

u H+1 =        B H+1     p 0 - H+1 i=0 B i B H+1 1-α m 0     +R H+1 CB H+1        (10)
Of course, if there is a limit N max < ∞, we take the minimum between the expression above and N max .

R H+1 is the unused leftover during the previous cycle in other words:

R H+1 = B H p 0 - H i=0 B i B H 1-α m 0 + R H -B H+1 C (11)
The big first term in the numerator of Eq. ( 10) is the total production minus the total maintenance. And at the denominator, we take into account the fact the B H+1 extremities share this total resource and C is the cost of a child. This formula only works as long as the tree grows "steadily", meaning that if along the flux up some branches could not pay maintenance or if the extremities do not have enough to create a child ( u H+1 = 0 ), then [START_REF] Bohn | [END_REF] and ( 11) stop being predictive. Conversely, as long as all the u H+1 give strictly positive integer, we can be sure no branches died along "flux up" (otherwise the numerator in (10) would be negative and so would u H+1 ). In other words, the moment u H+1 ≤ 0 is the moment the tree has stopped growing "steadily" and we note H f the last height/generation [START_REF] Bohn | [END_REF] is valid. H f may or may not be infinite but we are particularly interested in the infinite case.

6.2.1

The case R i = 0: the tree does not keep its leftover resource Now we put the term R H+1 to zero and since u H = B H+1 /B H , we obtain a simple expression to get the successive B H+1 .

B H+1 = B H •min       p 0 - H i=0 Bi B H 1-α m 0 C    , N max    (12)
To simplify things we introduce

W n = n i=0 (B i /B n ) 1-α . W n+1 = B n+1 B n+1 1-α + n i=0 B i B n+1 1-α = 1+ n i=0 B i B n+1 1-α .
Then we replace the B n+1 term at the denominator by the expression in (12) which replaces B n+1 by B n multiplied by a factor. This term B n allows us to recover a W n , and we end up with a "simple" recurrence.

f 0 (x)
x Fig. 6. Plot of f0, f1 and I for Nc max = 6 and α = 0.5. For 0 < α < 1, the f functions are monotonous increasing and f1 intersects twice the diagonal line. Here, a0 = 10 and b0 = 2.

W n+1 = 1 + W n 1 (min ( a 0 -W n • b 0 , N max )) 1-α
Where a 0 = p 0 /C and b 0 = m 0 /C with W 0 = 1 (13) Just as with Eq. ( 10), this expression is valid only until n = H f which corresponds to the first integer that yields a 0 -W H f • b 0 ≤ 0. Since we are interested in trees that can grow infinitely we limit our study to 0 < α < 1, and in this case H f may be infinite for some values of a 0 and b 0 . H f = ∞ is equivalent to W n (equivalently B n ) being defined for any natural number n. And if these two sequences are defined for every n we may want to look if they converge or not. First we need to establish when do we get an infinite sequence ( H f = ∞ ). We also remind a 0 , b 0 > 0. p 0 > m 0 is also assumed because otherwise our tree will stop growing at H f = 1.

We define the function f 0 such that Eq. ( 13) simply becomes a discrete dynamical system

W n+1 = f 0 (W n ) and f 1 (x) = 1 + x/ (a 0 -x • b 0 )
1-α , that we can study with classical tools. In fact f 1 is a simplified and a more easy-to-study version of f 0 . Their behavior with respect to I : x → x provides us with the behavior of (W n ). The figure (6) provides us with the behavior of these functions with respect to one another. Before moving to study this dynamical system, we point out a few things about f 0 , f 1 and W n :

-The domains of definition of our functions are D(f 1 ) = (-∞, a 0 /b 0 ) and D(f 0 ) = (-∞, (a 0 -1)/b 0 ]. The functions being positive and the initial condition being

W 0 = 1, only R + ∩ D(f 1 ) and R + ∩ D(f 0 ) interest us.
The sequence (W n ) n stops when f 0 (W n ) escapes the domain of definition of f 0 and the term H f is the last defined one. -Given a 0 , b 0 > 0 and 0 < α < 1, both f 0 and f 1 are strictly increasing in the domains we are interested in.

Then we can show that even for a non-continuous function like f 0 , we end up with an attractor which is a simple fixed point. (see Appendix 7 for details) this means that the sequence W n converge in most cases, it is indeed possible that f 0 is not below the line y = x. Now let us suppose we do have a fixed point. We call this fixed point/ limit W * . But

u n = min ([a 0 -W n b 0 ] , N max ) , (14)
so (u n ) n converges too. If we assume N max > a 0 -W * b 0 then u * = a 0 -W * b 0 . So if we fix a 0 and b 0 , the limits of both (W n ) and (u n ) n will be tied together. The conclusion, which is observable through the simulations, is that when we negate the reserve R = 0 then infinite trees give sequences (u i ) that ends with an infinite succession of u * .

In conclusion, with the dynamics and strategy we described (there is no external threat and the branches are acting brainlessly), infinite trees only appear for some values of the parameters and they end up, after a certain time, growing in a very regular way: the number of extremities is multiplied by u * each generation and we end up with a self-similar tree.

6.2.2 Case with the reserve on: back to equation [START_REF] Bohn | [END_REF].

First we can make the following remark: even with the reserve turned on, the results about the necessary conditions to get an infinitely growing tree described earlier should, at least to some extent, hold true. However with the addition of the reserve, we observe (from numerical calculations) that instead of having a sequence (u i ) i that ends with an infinite succession of u * , we may sometimes have a periodicity: the sequence will end up oscillating between u * 1 = a 0 -W * b 0 and u * 2 = a 0 -W * b 0 . There does not seem to be limits for the length of the periods we find. But, the period we get, vary progressively in function of the parameters a 0 and b 0 (Fig. 8). So while we have a less repetitive growth than in the case R = 0, the tree grows in a very regular manner and do not seem to change erratically its growth pattern when we vary slowly the parameter values.

Conclusion and perspectives

In this paper, we have presented a set of basic rules describing a group of bio-inspired dynamical systems focused on resource distribution and allocation. By completing these basic rules with specific evolution rules, we 3), while the one on the bottom-right corner encoded in white is the period (4). We can see that between (3) and ( 4), we go through a multitude of bands of different color. These periods are a succession of 3 and 4 but some have a bigger proportion of 4 or 3: the one near the (4) region has more 4 and vice versa. The large band in gray on the middle represents the period (3,[START_REF] Ford | Flows in Networks[END_REF] and it separates regions where the proportion of 4 is larger from the ones with more 3. There appears to be some selfsimilar patterns at the frontier of two "large" domains. There, we made a graph around a0 = 4 but if we make one around an other integer a0 = N , the main difference will be that the periods are composed of the integers N and N -1 but other than that we will see the same pattern of bands. construct a growing tree-like system we can simulate and study. However, even with only the basic rules and no evolution rules, we have shown it was possible to get a few results about the possible topologies of a tree-network that possesses sources at its extremities and have a maintenance cost for each node which increases the closer the node is to the root (Leonardo's rule). Then, we specified the evolution rules and studied the resulting dynamical system. By studying a very simple version of the system, we could established a few theoretical and numerical results that can be useful to lay some ground for future works. We emphasize as well that during our study we initially spend a lot of time scanning parameter space for the case α = 1, in this situation, as already discussed, most trees end up dying, however for some carefully chosen parameters we can end up with what would be an like an unstable fixed point, a tree that simply stops growing and lives forever. Nothing really interesting came out of this thorough numeric study which is why we did not present it in this paper, however these unstable tree allowed us to test the consequence of how the redistribution of resources among children affects the tree. The results are displayed in Fig. 9, we start from a perturbation of an unstable tree by adding three extra branches, this will lead to a depletion of the reserves and a death of the tree, but as illustrated we can see that depending on how the To be more precise the tree with rewarding strategy eternally oscillates between forms similar to the ones displayed in the top and middle figures, while the bottom tree dies shortly after, being unable to support its structure, this would lead to indicate a survival of the fittest strategy would be best in the case of this individual tree.

redistribution is made, the tree can be more resilient depending on the redistribution chosen . The formula that dictates the proportion a branch/parent will distribute its resource to each child depends on two elements: 1) the need of the child which is composed of its maintenance cost plus the cost of growing to reach the required volume (Leonardo's rule) 2) the amount of resource the child gave to the parent during "flux down". More specifically, if we note e i the need of a child i and c i its contribution during flux down then the formula is Ze r1 i c r2 i . With r 1 and r 2 parameters and Z the normalization constant j e r1 j c r2 j (sum over all the children j of this particular parent), and thus if r 1 > r 2 the tree redistribute resources based on needs, while the r 1 < r 2 is based more on reward. The results displayed in Fig. 9 clearly show that for the considered case, a strategy based on reward keeps the tree alive, while the one based on needs ends up killing the tree. This leaves us with the perspective of this work: we dubbed this paper part 1 kinematics, as no external forces or interaction between the branches besides redistribution of the resources was taken into account, so only the selfsustained kinematics of the tree were taken into account.

As a first perspective we want to embed this growth into real space, which will add occupation constraints on the new branches and some exclusion rules of available space to grow new offspring. This will lead to some interactions between the branches and we expect that the resulting dynamics will be greatly affected by this. Another aspect of future work, will be to influence the role of the redistribution parameter and its possible influence on the resilience of the tree structure as well as its overall shape like its extremities (foliage) when embedding it in real space. A comparison would be then possible with the already existing a attempts to obtain realistic looking tree shape from relatively simple rules [29,30]. Thus, with the embedding in space, we could examine whether our approach with resource distribution and allocation could yield similar results. Appendix A

Let us show that even for a non-continuous function like f 0 , we end up with an attractor which is a simple fixed point. Indeed, any real-valued function g strictly increasing defined on [x, y] such that g(x) > x and g(y) < y, g has a fixed point on [x, y]. Likewise with only these constraints on g , if we define u n+1 = g(u n ) with u 0 = x as initial condition then the sequence is define for any n and converges. We can apply this result to the function f 0 and the sequence (W n ) n with W 0 = 1. First, we may want to set a 0 > b 0 + 1 so that W 1 = f 0 (1) is defined. And if it is defined then it is obvious f 0 (1) > 1. Therefore, the last step is to prove the existence of a point z such that f 0 (z) < z, so we can deduce that (W n ) n is an infinite sequence and converges. On the other hand, if f 0 is always above the line y = x + then it is possible to show that (W n ) n will terminate at some integer H f < ∞. Since f 0 ≥ f 1 (cf. figures ( 6) and ( 7) ), studying f 1 can give us a sufficient condition for H f to be finite: having f 1 -I > would be that sufficient condition.

Let us study f 1 on D = R + ∩ D(f 1 ). Assuming again a 0 > b 0 + 1 > 1 and 0 < α < 1, we have both f 0,1 (1) > 1. Now we want to see when would f 1 -I > or not i.e. when we can not or can find y ∈ D such that f 1 (y) < y. The derivative of (f 1 -I) is: (15) Its sign is determined by its numerator that we will call P (x); P is defined on D(f 1 ).

(f 1 -I) (x) = ( 1
P (x) = b 0 (-α + (2 -α)(a 0 -b 0 x) 1-α ) ( 16
)
So P is, at first, strictly increasing until it reaches its maximum value at x 1 ≡ 1 b0 a 0 -α 2-α α then becomes strictly decreasing. In addition to that, because 2 -α > 1 we have P (-∞) < 0 and at the limit x → a 0 /b 0 we have P (x) > 0. So we can now establish a variation table for P and from there we deduce that there is a unique x 0 ∈ D(f 1 ) such that P (x 0 ) = 0 and P (x) < 0 for x < x 0 while P (x) > 0 for x > x 0 . We can even specify a bit more the value of x 0 : using the fact that a 0 > 1 we get P (0) = a 0 -a 2-α 0 < 0, furthermore P (x 1 ) > 0. The conclusion from these two facts is x 0 ∈ (0, x 1 ). (Remark: x 1 > a0-1 b0 because α 2-α α < 1 and a 0 > b 0 + 1 > 1.) From the sign of P , we finally deduce the variation table of f 1 -I. In conclusion, f 1 -I reaches its minimum at x 0 ∈ (0, x 1 ); x 0 = a0 b0 -X0 b0 where X 0 is the (unique) root in R + of a 0 (1 -α) + αX -X 2-α . If f 1 (x 0 ) -x 0 = with > 0 then H f is finite since it would imply f 1 -I ≥ as x 0 was the minimum of f 1 -I. So f 1 (x 0 ) -x 0 = is the sufficient (and probably nearly necessary) condition for H f < ∞ and it only requires a numerical determination of X 0 then a calculation of f 1 (x 0 ) = 1 + x 0 /X 0 .

Up to now, we saw that f 0 being below the line y = x for some real in D(f 0 ) implies the existence of a fixed point and H f = ∞. Then we have managed to find some sufficient condition for it not to be the case. Now let us suppose we do have a fixed point. Since f 0 (1) > 1 the fixed points should be located after 1. Considering how f 1 vary, it does not have more than 2 fixed points. We will assume the same for f 0 ; (W n ) n will then converge toward the smallest fixed point of f 0 that we will assume to be equal to the smallest fixed point of f 1 (both should be very close from each other).

Fig. 1 .

 1 Fig. 1. Acyclic connected networks represent trees. To illustrate the terminology we use in the paper: node T is the trunk, it has two children (node A and B which have two children).The tree has 5 extremities and is of height 3. And if the tree follows Leonardo's rule then the volume of T will be 5 units of volume, A will be 2 units and B will be 3 units.

Fig. 3 .

 3 Fig.3. A 3D rendering of a tree created from the model. The angle of the branches were chosen to make the drawing visually appealing but in our model trees exist as pure graphs free of spacial constraint.

Fig. 7 .

 7 Fig. 7. Using different parameters for a0 and b0 than in figure (6), we can elevate the graphs of f0,1 above the line y = x. In such case, we are sure that (Wn) is finite since one Wi will go beyond the domain of definition of f0,1.

Fig. 8 .

 8 Fig.8. The sequences (ui)i become periodic after a few terms.Here, for α = 0.5, we noted down the period we get for different values of a0 and b0 and we found 36 different periods. Each of them is encoded by one color. The period on the top-left corner encoded in black is the period composed only of 3, noted (3), while the one on the bottom-right corner encoded in white is the period (4). We can see that between (3) and (4), we go through a multitude of bands of different color. These periods are a succession of 3 and 4 but some have a bigger proportion of 4 or 3: the one near the (4) region has more 4 and vice versa. The large band in gray on the middle represents the period(3,[START_REF] Ford | Flows in Networks[END_REF] and it separates regions where the proportion of 4 is larger from the ones with more 3. There appears to be some selfsimilar patterns at the frontier of two "large" domains. There, we made a graph around a0 = 4 but if we make one around an other integer a0 = N , the main difference will be that the periods are composed of the integers N and N -1 but other than that we will see the same pattern of bands.

Fig. 9 .

 9 Fig. 9. A representation of three trees. The node at the bottom is the trunk. An eternal tree is perturbed adding 3 extremal branches (Top). Consequences of redistribution rules are shown in the middle and the bottom figure. The "rewarding" strategy ( r1 = 0.3 and r2 = 1) leads to an eternal tree (middle one) while the "altruist" strategy (r1 = 2 and r2 = 1) leads to a tree that ends up dying (bottom figure) . The values for the other parameters are: p0 = 400, m0 = 50, Cv = 10, C = 70 and α = 1.To be more precise the tree with rewarding strategy eternally oscillates between forms similar to the ones displayed in the top and middle figures, while the bottom tree dies shortly after, being unable to support its structure, this would lead to indicate a survival of the fittest strategy would be best in the case of this individual tree.

  Ackowledgements X.L. thanks Christophe Eloy, for fruitful discussions and contributions especially in the early development of the model, and its numerical implementation. O.B. and X.L. would like as well to thank Bruno Moulia, Eric Badel and André Lacointe for encouragements and useful suggestions. The project leading to this publication has received funding from Excellence Initiative of Aix-Marseille University -A*MIDEX, a French ?Investissements d?Avenir? programme. It has been carried out in the framework of the Labex MEC. We also acknowledge support from the CNRS (Mission pour l?interdisciplinarité, project AR-BRE)

  -α)b 0 x + a 0 -b 0 x -(a 0 -b 0 x) 2-α (a 0 -b 0 x) 2-α = a 0 (1 -α) + α(a 0 -b 0 x) -(a 0 -b 0 x) 2-α (a 0 -b 0 x) 2-α.