Yujie Wang

Yaozu Liu

Hui Li

Xinyu Guan

Ming Xue

Yushan Yan

Valentin Valtchev

Shilun Qiu

Qianrong Fang

Three-Dimensional Mesoporous Covalent Organic Frameworks through Steric Hindrance Engineering

The development of three-dimensional covalent organic frameworks (COFs) with large pores and high surface areas is of great importance for various applications. However, it remains a major challenge due to the unavoidable structural interpenetration and pore collapse after the removal of guest species sitiated in the pores. Herein, we report for the first time a series of 3D mesoporous COFs through a general strategy of enhanced steric hindrance. By using methoxy-modified monomer and increasing methyl groups of linkers, these 3D COFs can be obtained successfully as exclusively non-interpenetrated diamondoid structures, permanent mesopores (up to 26.5 Å), and high surface areas (> 3000 m 2 g -1 ), which are far superior to those of reported conventional COFs with the same topology. This work thus opens a way to create 3D large-porous COFs for potential applications in adsorption and separation of large inorganic, organic, and biological species.

Covalent organic frameworks (COFs) are fascinating crystalline porous polymers composed of light elements (typically H, B, C, N, and O) and linked by reversible covalent bonds. [START_REF] Côté | Crystalline, Covalent Organic Frameworks[END_REF][START_REF] Feng | Covalent Organic Frameworks[END_REF][START_REF] Ding | Covalent Organic Frameworks (COFs): from Design to Applications[END_REF][START_REF] Diercks | The Atom, the Molecule, and the Covalent Organic Framework[END_REF] They feature well-defined pores and high specific surface areas and thus hold great promises in a variety of potential uses, such as adsorption and separation, [START_REF] Han | III Covalent Organic Frameworks as Exceptional Hydrogen Storage Materials[END_REF][START_REF] Kuhn | Covalent TriazineBased Frameworks Prepared by Ionothermal Synthesis[END_REF][START_REF] Wang | Exfoliation of Covalent Organic Frameworks into Few-layer Redox-Active Nanosheets as Cathode Materials for Lithium-Ion Batteries[END_REF] organic electronics, [START_REF] Wan | A Belt-Shaped, Blue Luminescent, and Semiconducting Covalent Organic Framework[END_REF][START_REF] Dogru | A Photoconductive Thienothiophene-Based Covalent Organic Framework Showing Charge Transfer Towards Included Fullerene[END_REF][START_REF] Bertrand | Thiophene-Based Covalent Organic Frameworks[END_REF][START_REF] Du | Ionic Covalent Organic Frameworks with Spiroborate Linkage[END_REF][START_REF] Li | Three-Dimensional Tetrathiafulvalene-Based Covalent Organic Frameworks for Tunable Electrical Conductivity[END_REF] heterogeneous catalysis, [START_REF] Ding | Construction of Covalent Organic Framework for Catalysis: Pd/COF-LZU1 in Suzuki-Miyaura Coupling Reaction[END_REF][START_REF] Peng | Synthesis of a Sulfonated Two-Dimensional Covalent Organic Framework as an Efficient Solid Acid Catalyst for Biobased Chemical Conversion[END_REF][START_REF] Vyas | A Tunable Azine Covalent Organic Framework Platform for Visible Light-Induced Hydrogen Generation[END_REF][START_REF] Sun | Flexibility Matters: Cooperative Active Sites in Covalent Organic Framework and Threaded Ionic Polymer[END_REF][START_REF] Wang | Homochiral 2D Porous Covalent Organic Frameworks for Heterogeneous Asymmetric Catalysis[END_REF] and some others. [START_REF] Oh | A Cryogenically Flexible Covalent Organic Framework for Efficient Hydrogen Isotope Separation by Quantum Sieving[END_REF][START_REF] Zhou | One-Step Construction of Two Different Kinds of Pores in a 2D Covalent Organic Framework[END_REF][START_REF] Chandra | Phosphoric Acid Loaded Azo (-N=N-) Based Covalent Organic Framework for Proton Conduction[END_REF][START_REF] Jin | Two-Dimensional sp2 Carbon-Conjugated Covalent Organic Frameworks[END_REF][START_REF] Guan | Chemically Stable Polyarylether-Based Covalent Organic Frameworks[END_REF][START_REF] Guo | Fast Ion Transport Pathway Provided by Polyethylene Glycol Confined in Covalent Organic Frameworks[END_REF] Over the past decade, most of COF materials, however, were exploited based on two-dimensional (2D) structures with characteristic eclipsed AA stacking. Three-dimensional (3D) COFs have recently attracted widespread attentions owing to their unique pores and excellent properties. [START_REF] El-Kaderi | Designed Synthesis of 3D Covalent Organic Frameworks[END_REF][START_REF] Uribe-Romo | A Crystalline Imine-Linked 3D Porous Covalent Organic Framework[END_REF][START_REF] Beaudoin | Constructing Monocrystalline Covalent Organic Networks by Polymerization[END_REF][START_REF] Liu | Weaving of Organic Threads into a Crystalline Covalent Organic Framework[END_REF][START_REF] Li | Three-Dimensional Covalent Organic Frameworks with Dual Linkages for Bifunctional Cascade Catalysis[END_REF][START_REF] Li | ThreeDimensional Ionic Covalent Organic Frameworks for Rapid, Reversible, and Selective Ion Exchange[END_REF][START_REF] Ma | Single-Crystal X-Ray Diffraction Structures of Covalent Organic Frameworks[END_REF][START_REF] Lu | Postsynthetic Functionalization of Three-Dimensional Covalent Organic Frameworks for Selective Extraction of Lanthanide Ions[END_REF][START_REF] Yan | Three-dimensional Salphen-based Covalent-Organic Frameworks as Catalytic Antioxidants[END_REF] Notably, the urgent issues emerge for the construction of 3D large-porous COFs, including inevitable structural interpenetration for restricting the pore size and pore shrinkage or collapse after the guest removal. For example, Wang and Zhang et al. have proven the reversible crystal-transformation and dynamic response in an interpenetrating 3D COF, LZU-301, which has pore size stretchable from 5.8 × 10.4 Å to 9.6 × 10.4 Å by the induction of guests. [START_REF] Ma | A Dynamic Three-Dimensional Covalent Organic Framework[END_REF] Recently, by extending the linear linkers, we obtained the most interpenetrated framework for 3D COFs reported, 3D-IL-COF-3 with 11-fold interpenetrated net, exhibiting a small pore size of 12.4 Å and relatively low specific surface area of 870 m 2 /g. [START_REF] Guan | Ambient Temperature and Pressure Ionothermal Synthesis of Three-Dimensional Covalent Organic Frameworks[END_REF] Thus, the expansion of 3D COFs with large pores and high surface areas is of primary importance but remains a tremendouschallenge.

Herein, we report a series of 3D mesoporous COFs (3D-Meso-COFs) with unusual non-interpenetrated diamondoid (dia) topology based on the inovative enhanced steric hindrance strategy. By employing methoxy-modified monomer and changing substituents (methyl groups) of linkers, these 3D-Meso-COFs can be turned successfully from pore shrinkage to partial and full pore maintenance. Remarkably, one of the obtained COF material shows permanent mesopore with a size of 26.5 Å and impressive specific surface area of 3023 m 2 g -1 , which are much higher than that of the convention interpenetrating COF with the same network. To the best of our knowledge, this study Scheme 1. Strategy for preparing 3D mesoporous COFs a a (a) The interpenetrated microporous framework constructed from undecorated edges. (b) The non-interpenetrated mesoporous framework built up from partially decorated edges, in which the pore shrinkage or collapse will happen after the guest removal. (c) The non-interpenetrated mesoporous framework created by highly decorated edges, in which the pores can be maintained entirely after the guest removal.

represents the first case of 3D mesoporous COF with non-interpenetrated dia topology.

Our strategy for the preparation of 3D-Meso-COFs is based on the steric hindrance by employing various linkers. As shown in Scheme 1three possible situations can occur: (a) the interpenetrated microporous framework constructed from undecorated edges; (b) the non-interpenetrated mesoporous framework built up from partially decorated edges, in which the pore shrinkage or collapse will happen after the guest removal; and (c) the non-interpenetrated mesoporous framework created by highly decorated edges, in which the pores can be maintained entirely after the removal of guests. To implement this strategy, we firstly designed a novel methoxy-modified monomer, 2,2',7,7'-tetramethoxy-9,9'-spirobi[fluorene]-3,3',6,6'-tetracarbald ehyde (TMSFTA), as a tetrahedral building unit (Figure 1a). Further, we chose three linear linkers, 4,4'-diaminobiphenyl (DABP, Figure 1b) without any substituents, 4,4'-diamino-2,2'-dimethylbiphenyl (DADMBP, Figure 1c) with two methyl groups, and 4,4'-diaminp-2,2',6,6'-tetramethylbiphenyl (DATMBP, Figure 1d) with four methyl groups. The condensation of TMSFTA and DABP, DADMBP or DATMBP produced a series of 3D COFs, termed JUC-550, JUC-551 and JUC-552 respectively (JUC = Jilin University China, Figure 1e-g). In light of the linking of tetrahedral and linear building units and their steric hindrance, these structures are expected to be the non-interpenetrated dia nets (Figure 1h and1i). 35 Typically, the syntheses were carried out by suspending TMSFTA with DABP, DADMBP or DATMBP into the solvent of dioxane or mesitylene in the presence of acetic acid, followed by heating at 120 °C for 3 days to give crystalline solids at high yields of 85% for JUC-550, 83% for JUC-551, and 87% for JUC-552, respectively. These products were stable in various organic solvents, such as acetone, hexanes, tetrahydrofuran (THF), N,N-dimethylformamide (DMF), N-methyl-2-pyrrolidone (NMP), and dimethyl sulfoxide (DMSO).

The structural characteristics of new COFs were determined by the complementary method. The crystal morphology was observed by scanning electron microscopy (SEM), which shows isometric crystals (Figures S1-3). Fourier transform infrared (FT-IR) spectra exhibited new peaks around 1617 cm -1 for JUC-550, 1612 cm -1 for JUC-551, and 1617 cm -1 for JUC-552, corresponding to typical C=N stretching vibrations (Figures S4-6). Furthermore, the solid state [START_REF] Ding | Construction of Covalent Organic Framework for Catalysis: Pd/COF-LZU1 in Suzuki-Miyaura Coupling Reaction[END_REF] C cross-polarization magic-angle-spinning (CP/MAS) NMR spectra confirmed the presence of carbon from the C=N bond at 158 ppm for JUC-550, 161 ppm for JUC-551, and 159 ppm for JUC-552, respectively (Figures S7-9). According to the thermogravimetric analysis (TGA), these COFs were thermally stable up to 350 °C under nitrogen (Figures S10-12). 224), (314) and (402) Bragg peaks of the same space group, respectively. The refinement results yielded unit cell parameters nearly equivalent to the predictions with good agreement factors (Rp = 2.58% and ωRp = 3.77% for JUC-550; Rp = 3.98% and ωRp = 5.06% for JUC-551; and Rp = 2.59% and ωRp = 3.81% for JUC-552). Some alternative structures, such as 2-and 3-fold interpenetrated dia nets, were also set up; however, their simulated PXRDs did not match to the experimental ones (Figures S13-24). Thus, the obtained COFs were proposed to have the expected mesoporous frameworks with non-interpenetrated dia topology (Figure 3).

It should be noted that some peaks after 2 ˃ 15° in JUC-550 are intense, while those of JUC-551 showed low intensity. Most important, the corresponding peaks of JUC-552 were very week revealing the impact of the linker on the crystalline framework. Furthermore, after activation under vacuum at 100 °C for 12 h, the first peak (2 = 3.35°) was disappeared in JUC-550, while showed a slight variation in JUC-551. However, JUC-552 exhibited enhanced intensity for its first peak (2 = 3.34°), which means its mesopores have been well retained. Notably, this process can be repeated after immersing these samples into organic solvents, such as acetone (Figures S25-27). These results demonstrated clearly that a guest-induced reversible crystal-structure transformation took place in JUC-550; [START_REF] Ma | A Dynamic Three-Dimensional Covalent Organic Framework[END_REF] however, the linker with high substituents, DATMBP, gives rise to JUC-552 with increased structural rigidity due to the strong steric hindrance, and thus can effectively eliminate the dynamic behavior of frameworks. 37 The influence of steric hindrance in activated 3D-Meso-COFs was further proved by nitrogen (N2) adsorption-desorption measurment at 77 K. As can be seen in Figure 4a, a rapid uptake at low relative pressure characteristic of micropore material was observed for JUC-550, which reveals a pore shrinkage. JUC-551 showed a rapid uptake at a low pressure of P/P0 < 0.05, followed by two other steps at P/P0 = 0.1 and 0.4, which means its partial pore shrinkage and the presence of different pores. On the contrary, JUC-552 had a high uptake at a low pressure of P/P0 < 0.05, followed by well define second step at P/P0 = 0.3, which is a typical characteristic of mesoporous materials. All curves of 3D-Meso-COFs showed a desorption hysteresis, which can be attributed to the effect of substituents and textural mesopores from the agglomeration of COF microcrystals. [START_REF] Fang | 3D Porous Crystalline Polyimide Covalent Organic Frameworks for Drug Delivery[END_REF] The Brunauer-Emmett-Teller (BET) equation was employed to evaluate the specific surface area of the series of COF materials. A strong correlation between the linker employed and the SBET was found, thus JUC-550, JUC-551 and JUC-552 showed 846 m 2 g -1 1728 m 2 g -1 , and 3023 m 2 g -1 , respectively. Furthermore, their pore-size distributions calculated by nonlocal density functional theory (NLDFT) demonstrated the micropores with the size of 10.6 Å for JUC-550 and two kinds of pores with the sizes of 14.0 Å and 22.0 Å for JUC-551. Notably, JUC-552 showed the mesopores with the dominant size of 26.5 Å, which is in good agreement with that of the proposed model (26.0 Å), and thus indicated that the pores can be maintained after the guest removal. As a comparison, we summarized the BET surface areas and pore sizes of current reported COFs with dia topology in Figure 5 and Table S1. Remarkably, the pore size of JUC-552 (26.5 Å) is far superior to those of reported materials, such as 4.0 Å for COF-300-C7, [START_REF] Ma | Observation of Interpenetration Isomerism in Covalent Organic Frameworks[END_REF] 12.0 Å for COF-320, [START_REF] Zhang | Single-Crystal Structure of a Covalent Organic Framework[END_REF] 13.0 Å for JUC-508, [START_REF] Yan | Three-dimensional Salphen-based Covalent-Organic Frameworks as Catalytic Antioxidants[END_REF] and 13.5 Å for 3D-OH-COF. [START_REF] Lu | Postsynthetic Functionalization of Three-Dimensional Covalent Organic Frameworks for Selective Extraction of Lanthanide Ions[END_REF] The BET surface area of JUC-552 (3023 m 2 g -1 ) is also much higher than those of reported conventional COFs, such as 517 m 2 g -1 for 3D-IL-COF-1, 34 864 m 2 g -1 for 3D-HNU-5, 41 1762 m 2 g -1 for COF-DL229, [START_REF] Wang | A 3D Covalent Organic Framework with Exceptionally High Iodine Capture Capability[END_REF] and 2403 m 2 g -1 for PI-COF-4. [START_REF] Fang | 3D Porous Crystalline Polyimide Covalent Organic Frameworks for Drug Delivery[END_REF] Undoubtly, these data revealed that the enhanced steric hindrance is a promising strategy for the construction of COF materials with large pores and high specific surface areas.

To further demonstrate the mesopores in 3D-Meso-COFs, we explored the inclusion of a large dye molecule, Rh6G, with a size of 16.0 Å. Typically, the samples of 3D-Meso-COFs were immersed in a solution of Rh6G (0.1 mM). The amount of Rh6G in the supernatant was monitored through UV-Vis spectrophotometry, and the characteristic absorbance at 526 nm was measured over a period of 60 mins. As showed in Figures S36 andS37, the results displayed a continuous decrease in the amount of Rh6G in solvated 3D-Meso-COFs due to their maintained mesopores. As a control, a similar experiment was carried out with COF-320, whose pore size (12.0 Å) is smaller than Rh6G. the experiment confirmed the absence of adsorption Rh6G on COF-320, no obvious reduction in absorbance was found (Figures S36 andS37).

In conclusion, we design and synthesize a series of 3D mesoporous COFs by a simple and effective strategy of enhanced steric hindrance. By utilizing methoxy-modified linker, 3D-Meso-COFs can be prepared as non-interpenetrated dia structures. The new strategy allows fine tuning the properties of COF mesoporosity by controlling the amount of substutient in the framework. Thus a new material, JUC-552 with stable mesopores (~ 26.5 Å) and outstanding specific surface areas (up to 3023 m 2 g -1 ), which are much higher than those of reported COFs with the same topology, was obtained. This research starts a new avenue to construct 3D large-porous COFs for various applications, such as adsorption and separation of biological, inorganic, and organic molecules with large sizes.

ASSOCIATED CONTENT Supporting Information

Synthetic procedures, SEM, FTIR, solid state [START_REF] Ding | Construction of Covalent Organic Framework for Catalysis: Pd/COF-LZU1 in Suzuki-Miyaura Coupling Reaction[END_REF] C NMR, TGA, BET plots, and pore size distribution. This material is available free of charge via the internet at http://pubs.acs.org.

Figure 1 .

 1 Figure 1. (a-d) Molecular structures of TMSFTA (a) as a tetrahedral building unit and DABP (b), DADMBP (c) and DATMBP (d) as a linear building unit, respectively. (e-g) 3D mesoporous COFs, denoted as JUC-550 to -552, constructed by the condensation reaction of TMSFTA and DABP, DADMBP or DATMBP. (h and i) Single dia network (h) and 3D non-interpenetrated dia structure (i) with decorated edges in these obtained COFs.

Figure 2 .

 2 Figure 2. PXRD patterns of JUC-550 (a), JUC-551 (b) and JUC-552 (c).

Figure 3 .

 3 Figure 3. Structural representations of JUC-550 (a), JUC-551 (b), and JUC-552 (c). C, blue; H, gray; N, green; O, pink. The spheres with different colors represent the pores in 3D-Meso-COFs. The structures of 3D-Meso-COFs were resolved by powder X-ray diffraction (PXRD) combined with structural simulations (Figure 2). After a geometrical energy minimization by using the Materials Studio software package 36 based on the non-interpenetrated dia net, the unit cell parameters were obtained (a = b = 27.5727 Å, c = 74.5020 Å and α = β = γ = 90° for JUC-550; a = b = 27.8446 Å, c = 74.6490 Å and α = β = γ = 90° for JUC-551; and a = b = 27.8435 Å, c = 74.6423 Å and α = β = γ = 90° for JUC-552). The simulated PXRD patterns were in good agreement with the experimental results. Full profile pattern matching (Pawley) refinements were carried out on the experimental data from solvated samples. Peaks at 3.35, 4.98, 6.74, 7.53, 10.07, 11.96 and 13.65° for JUC-550 correspond to the (101), (112), (200), (202), (301), (321) and (404) Bragg peaks of space group I41/a (No. 88); peaks at 3.38, 4.79, 5.08, 6.37, 6.82, 8.01, 9.71, 10.29, 11.23 and 13.07° for JUC-551 correspond to the (101), (103), (112), (200), (202), (213), (301), (224), (314) and (402) Bragg peaks of the same space group; and peaks at 3.34, 4.78, 5.07, 6.35, 6.79, 7.99, 9.68, 10.23, 11.29 and 13.01° for JUC-552 correspond to the (101), (103), (112), (200), (202), (213), (301), (224), (314) and (402) Bragg peaks of the same space group, respectively. The refinement results yielded unit cell parameters nearly equivalent to the predictions with good agreement factors (Rp = 2.58% and ωRp = 3.77% for JUC-550; Rp = 3.98% and ωRp = 5.06% for JUC-551; and Rp = 2.59% and ωRp = 3.81% for JUC-552). Some alternative structures, such as 2-and 3-fold interpenetrated dia nets, were also set up; however, their simulated PXRDs did not match to the experimental ones (FiguresS13-24). Thus, the obtained COFs were proposed to have

Figure 4 .

 4 Figure 4. N2 adsorption-desorption isotherms for JUC-550 (a), JUC-551 (b) and JUC-552 (c).

Figure 5 .

 5 Figure 5. A plot of BET surface area against pore size in reported COFs with dia topology.

ACKNOWLEDGMENT This work was supported by National Natural Science Foundation of China (21571079, 21621001, 21390394, 21571076, and 21571078), 111 project (B07016 and B17020), and the program for JLU Science and Technology Innovative Research Team. V.V. and Q.F. acknowledge the Thousand Talents program (China).

AUTHOR INFORMATION

Corresponding Author *qrfang@jlu.edu.cn Author Contributions § Y.W. and Y.L. contributed equally to this work.

Notes

The authors declare no competing financial interests.


Vincenzo Catrambone

Herwig Wendt

Riccardo Barbieri

Patrice Abry

Gaetano Valenza

Quantifying Functional Links between Brain and Heartbeat Dynamics in the Multifractal Domain: a Preliminary Analysis

Quantification of brain-heart interplay (BHI) has mainly been performed in the time and frequency domains. However, such functional interactions are likely to involve nonlinear dynamics associated with the two systems. To this extent, in this preliminary study we investigate the functional coupling between multifractal properties of Electroencephalography (EEG) and Heart Rate Variability (HRV) series using a channel-and time scale-wise maximal information coefficient analysis. Experimental results were gathered from 24 healthy volunteers undergoing a resting state and a cold-pressure test, and suggest that significant changes between the two experimental conditions might be associated with nonlinear quantifiers of the multifractal spectrum. Particularly, major brain-heart functional coupling was associated with the secondorder cumulant of the multifractal spectrum. We conclude that a functional nonlinear relationship between brain-and heartbeat-related multifractal sprectra exist, with higher values associated with the resting state.

Clinical relevance: The presented methodology could provide biomarkers for supporting the diagnosis of neurological, psychiatric, and cardiovascular disorders.

I. INTRODUCTION

The central nervous system and the cardiovascular system continuously interact and communicate through anatomical, biochemical, and functional links, and the resulting functional network has been referred to as the central-autonomic network (CAN) [START_REF] Valenza | The central autonomic network at rest: Uncovering functional mri correlates of time-varying autonomic outflow[END_REF]- [START_REF] Valenza | Uncovering complex central autonomic networks at rest: a functional magnetic resonance imaging study on complex cardiovascular oscillations[END_REF].

Focusing on Electroencephalography (EEG) and Heart Rate Variability (HRV) series, previous studies investigated CAN-related changes between different populations and experimental conditions, and model-based and model-free approaches have been proposed. While a model-based approach exploits a-priori knowledge or assumptions to try to describe the functional brain-heart information exchange, a modelfree approach focuses on synchronized brain and heartbeat events. Exemplarily, we recently proposed a model-based approach based on synthetic data generation for a functional directional brain-heart interplay assessment, which was evaluated with data gathered from healthy subjects undergoing relevant sympathovagal changes induced by a cold-pressor test [START_REF] Catrambone | Time-resolved directional brain-heart interplay measurement through synthetic data generation models[END_REF]. On the other hand, exemplary model-free approaches include the so-called Heartbeat-Evoked Potential (HEP) [START_REF] Schandry | From the heart to the brain: a study of heartbeat contingent scalp potentials[END_REF], the application of information theory methods applied to brain and heartbeat dynamics [START_REF] Faes | Information dynamics of brain-heart physiological networks during sleep[END_REF], and joint symbolic analysis [START_REF] Schulz | Central-and autonomic nervous system coupling in schizophrenia[END_REF].

The cardiovascular system shows a nonlinear and nonstationary behaviour due to several interactions with the respiratory and nervous systems [START_REF] Sunagawa | Dynamic nonlinear vago-sympathetic interaction in regulating heart rate[END_REF]. Consequently, heartbeat dynamics has been successfully modelled by exploiting the so-called multifractal (MF) framework analysis through the definition of a multifractal spectrum D(h) characterizing transient and local non-Gaussian structures that a single Hurst exponent, H, may not describe [START_REF] Ivanov | Multifractality in human heartbeat dynamics[END_REF]- [START_REF] Wendt | Wavelet p-leader non gaussian multiscale expansions for heart rate variability analysis in congestive heart failure patients[END_REF]. Similarly, EEG series show multiple local singular behaviours going beyond self-similarity and thus need a model comprising a collection of exponents H = h(t) to fully characterize its dynamics [START_REF] Vicsek | Fractal growth phenomena[END_REF], [START_REF] Takayasu | Fractals in the physical sciences[END_REF]. Exemplarily, MF properties were investigated in several EEG studies on, e.g., sleep [START_REF] Song | Multifractal analysis of sleep EEG dynamics in humans[END_REF], motor imagery [START_REF] Popivanov | Multifractality of decomposed EEG during imaginary and real visual-motor tracking[END_REF], and epileptic seizures [START_REF] Zhang | Multifractal analysis and relevance vector machinebased automatic seizure detection in intracranial EEG[END_REF], learning [START_REF] Rocca | Self-similarity and multifractality in human brain activity: a wavelet-based analysis of scale-free brain dynamics[END_REF], and functional brain connectivity [START_REF] Ciuciu | Interplay between functional connectivity and scale-free dynamics in intrinsic fmri networks[END_REF].

From an estimation viewpoint, the wavelet transform modulus maxima method [START_REF] Muzy | Multifractal formalism for fractal signals: The structure-function approach versus the wavelet-transform modulusmaxima method[END_REF], and the MF detrended fluctuation analysis [START_REF] Kantelhardt | Multifractal detrended fluctuation analysis of nonstationary time series[END_REF] were proposed to estimate D(h). Recently, we used a wavelet leader multifractal formalisms [START_REF] Wendt | Bootstrap for empirical multifractal analysis[END_REF] and its generalization using p-leaders [START_REF] Jaffard | p-exponent and p-leaders, part i: Negative pointwise regularity[END_REF], [START_REF] Leonarduzzi | p-exponent and p-leaders, part ii: Multifractal analysis. relations to detrended fluctuation analysis[END_REF] with non-gaussian expansion to characterize brain and heartbeat dynamics separately during a Cold-Pressor Test (CPT) [START_REF] Catrambone | Heartbeat dynamics analysis under coldpressure test using wavelet p-leader non-gaussian multiscale expansions[END_REF], [START_REF] Catrambone | Wavelet p-leader non-gaussian multiscale expansions for eeg series: an exploratory study on cold-pressor test[END_REF]. CPT comprises a strong thermal elicitation with consequent sympathovagal changes and was associated with significantly lower MF features derived from both EEG [START_REF] Catrambone | Wavelet p-leader non-gaussian multiscale expansions for eeg series: an exploratory study on cold-pressor test[END_REF] and heartbeat dynamics [START_REF] Catrambone | Heartbeat dynamics analysis under coldpressure test using wavelet p-leader non-gaussian multiscale expansions[END_REF] than resting state. Furthermore, the crucial role of preprocessing on heartbeat dynamics based on inhomogeneous point-process models has been demonstrated [START_REF] Catrambone | Heartbeat dynamics analysis under coldpressure test using wavelet p-leader non-gaussian multiscale expansions[END_REF], [START_REF] Valenza | Mortality prediction in severe congestive heart failure patients with multifractal point-process modeling of heartbeat dynamics[END_REF].

Despite the significance of the aforementioned evidences, the functional link between brain and heartbeat dynamics at a MF level has not been investigated yet. To this end, in this preliminary study we report on a methodological attempt aimed to combine brain-and heartbeat-derived multifractal features at different time scales using a channel-wise Maximal Information Coefficient (MIC) analysis. Experimental results were gathered from healthy volunteers undergoing resting state and CPT, and are shown in terms of p-values topographic maps.

II. MATERIALS AND METHODS

A. Experimental Setup

A 128-channel EEG and one-lead ECG were collected from thirty healthy volunteers (gender balanced, average age of 26.7) with a 500Hz sampling rate. Signals were gathered through a Geodesic EEG Systems 300 (Electrical Geodesics, Inc.). All subjects declared to be healthy and right-handed, and signed an informed consent. The experimental procedure was approved by the local ethical committee. Data with significant artifacts gathered from six participants were not considered for further analyses.

Subjects underwent a 3 min resting state session followed by a 3 min Cold-Pressor Test (CPT) session with their left hand in iced water, and were free to stop the session in case of pain. Further information can be found in [START_REF] Catrambone | Time-resolved directional brain-heart interplay measurement through synthetic data generation models[END_REF].

B. EEG and ECG Preprocessing

EEG preprocessing was applied to a subset of 90 channels over the scalp and followed the so-called HAPPE procedure, whose details are extensively reported in [START_REF] Gabard-Durnam | The harvard automated processing pipeline for electroencephalography (happe): standardized processing software for developmental and high-artifact data[END_REF]. For each EEG channel, the power spectral density (PSD) was then derived within five standard frequency bands: [START_REF] Wendt | Wavelet p-leader non gaussian multiscale expansions for heart rate variability analysis in congestive heart failure patients[END_REF][START_REF] Vicsek | Fractal growth phenomena[END_REF][START_REF] Takayasu | Fractals in the physical sciences[END_REF][START_REF] Song | Multifractal analysis of sleep EEG dynamics in humans[END_REF][START_REF] Popivanov | Multifractality of decomposed EEG during imaginary and real visual-motor tracking[END_REF][START_REF] Zhang | Multifractal analysis and relevance vector machinebased automatic seizure detection in intracranial EEG[END_REF][START_REF] Rocca | Self-similarity and multifractality in human brain activity: a wavelet-based analysis of scale-free brain dynamics[END_REF][START_REF] Ciuciu | Interplay between functional connectivity and scale-free dynamics in intrinsic fmri networks[END_REF][START_REF] Muzy | Multifractal formalism for fractal signals: The structure-function approach versus the wavelet-transform modulusmaxima method[END_REF][START_REF] Kantelhardt | Multifractal detrended fluctuation analysis of nonstationary time series[END_REF][START_REF] Wendt | Bootstrap for empirical multifractal analysis[END_REF][START_REF] Jaffard | p-exponent and p-leaders, part i: Negative pointwise regularity[END_REF][START_REF] Leonarduzzi | p-exponent and p-leaders, part ii: Multifractal analysis. relations to detrended fluctuation analysis[END_REF][START_REF] Catrambone | Heartbeat dynamics analysis under coldpressure test using wavelet p-leader non-gaussian multiscale expansions[END_REF][START_REF] Catrambone | Wavelet p-leader non-gaussian multiscale expansions for eeg series: an exploratory study on cold-pressor test[END_REF][START_REF] Valenza | Mortality prediction in severe congestive heart failure patients with multifractal point-process modeling of heartbeat dynamics[END_REF][START_REF] Gabard-Durnam | The harvard automated processing pipeline for electroencephalography (happe): standardized processing software for developmental and high-artifact data[END_REF][START_REF] Citi | A real-time automated point-process method for the detection and correction of erroneous and ectopic heartbeats[END_REF][START_REF] Mallat | A wavelet tour of signal processing[END_REF] and γ ∈ [30-70] (all expressed in Hz). The PSD was estimated by exploiting the well-known Welch method with a Hamming window of 2 sec with a moving step of 0.25 sec, thus obtaining a 4 Hz sampling rate for the time-frequency representation.

δ ∈ [1-4], θ ∈ [4-8], α ∈ [8-12], β ∈ [
The ECG was processed to derive heartbeat dynamics. Rwave events were automatically identified through the wellknown Pan-Tompkins algorithm, and eventual algorithmic and physiological artefacts were corrected using a previously developed method based on point-process statistics [START_REF] Citi | A real-time automated point-process method for the detection and correction of erroneous and ectopic heartbeats[END_REF]. According to our previous evidence [START_REF] Valenza | Mortality prediction in severe congestive heart failure patients with multifractal point-process modeling of heartbeat dynamics[END_REF], instantaneous firstorder statistics for heartbeat dynamics were retrieved from the instantaneous mean of an inverse-Gaussian probability density function predicting R-wave events in time, µ RR (t), in line with an inhomogeneous point-process framework [START_REF] Valenza | Mortality prediction in severe congestive heart failure patients with multifractal point-process modeling of heartbeat dynamics[END_REF].

C. Multifractal analysis

The wavelet spectrum of self-similar models is defined as:

S d X (j, q = 2) = 1 n j nj k=1 |d X (j, k)| q K2 jqH (1)
with a power-law exponent controlled by the well-known Hurst parameter H. Here, n j is the number of d X (j, k) available at scale 2 j , and d X (j, k) = ψ j,k |X is a pattern of oscillation reference narrow in the time and frequency domains with mother wavelet ψ [START_REF] Mallat | A wavelet tour of signal processing[END_REF].

The wavelet spectrum is bounded to the Fourier spectrum, and H depends on linear data properties across frequency bands [START_REF] Doret | Multifractal analysis of fetal heart rate variability in fetuses with and without severe acidosis during labor[END_REF]. On the other hand, multifractal models build upon self-similarity by replacing the static parameter H with a collection of local exponents H = h(t), and the multifractal spectrum D(h) is then defined as the temporal repartition of h(t) [START_REF] Doret | Multifractal analysis of fetal heart rate variability in fetuses with and without severe acidosis during labor[END_REF], [START_REF] Wendt | Bootstrap for empirical multifractal analysis[END_REF]. The estimation of D(h) from real data requires to impose q = 2 in (1) with positive and negative interval of moments q, and the use of wavelet pleaders X are defined as p norms of wavelet coefficients in a narrow temporal neighbourhood λ j,k = [k2 j , (k +1)2 j ), with 3λ j,k = m{-1,0,1} λ j,k+m , over all finer scales. Estimates from a multifractal spectrum D(h) can be derived from the log-leaders cumulants as follows:

(p) X : (p) X (j, k) = 2 j λ ⊂3λ j,k 2 -j |d X (λ )| p 1/p (2)
C (p) m (j) ≡ Cum m log (p) X (j) c 0 m + c m log(2 j ) (3)
More specifically, the first-order cumulant c 1 identifies the H-related linear properties and the D(h) mode, while the higher-order order cumulants c 2 , c 3 , c 4 represent nonlinear features from the multifractal spectrum such as width, asymmetry, and kurtosis. Further details are reported in [START_REF] Wendt | Bootstrap for empirical multifractal analysis[END_REF]. In this preliminary study we limited our analysis to the firstand second-order cumulants c 1 and c 2 .

D. Maximal Information Coefficient

MIC is an estimator of the linear and nonlinear coupling between pair of samples and is directly defined from the scatterplot of the two paired series, x and y [START_REF] Reshef | Detecting novel associations in large data sets[END_REF]. A grid with any number of rows and columns can be superimposed to the scatterplot. Considering G nx×ny as the set of all the possible partitions of that grid with n x number of rows and n y number of columns, the algorithm calculates the associated mutual information to all elements g ∈ G nx×ny and finds the normalized maximum among G as follows:

m nx×ny = max {I g } g∈Gn x×ny log min{n x , n y } (4)
MIC is then calculated as the maximal m nx×ny over all the possible pairs (n x , n y ). It has been found that MIC may also be estimated as: M IC(x, y) = max nxny<B {m nx×ny }, with B empirically defined as B = n 0.6 . A fully detailed description is reported in [START_REF] Reshef | Detecting novel associations in large data sets[END_REF].

E. Quantification of Functional Brain-Heart Interplay

For each of the twenty-four subjects and for each of the two experimental sessions, the two MF cumulants c 1 and c 2 were estimated separately from brain and heartbeat Fig. 2: Topographic distribution of MIC differences between CPT and resting phase for the first-order cumulant of the MF spectrum c 1 calculated between five EEG frequency powers and the first-order moment µ RR derived from point-process models for heartbeat dynamics. Fig. 3: Topographic distribution of MIC differences between CPT and resting phase for the second-order cumulant of the MF spectrum c 2 calculated between five EEG frequency powers and the first-order moment µ RR derived from point-process models for heartbeat dynamics. dynamics. Specifically, from the brain side, we calculated c 1 and c 2 for each of the 90 EEG channels and the 5 PSD time series at the standard frequency bands. From the heart side, we derived c 1 and c 2 associated with the µ RR series. For each MF cumulants (c 1 and c 2 ) separately, a MIC was calculated group-wise between vectors whose length was of 216 elements, resulting from the number of subjects (i.e. 24) times the number of time scales (i.e. 9). Thus, a functional brain-heart interplay was quantified as M IC c1 (x, y) and M IC c2 (x, y), where x represents estimates from the heart side (e.g., c 1 (µ RR )) and y represents for estimates from the brain side at a specific frequency band and EEG channel (e.g., c 1 (P SD α,F p2 )). A schematic description of the analysis pipeline is reported in Figure 1.

Finally, separately for c 1 and c 2 and for each frequency band, statistically significant differences between the topographic distributions of MIC values (across the 90 EEG channels) extracted for the two experimental conditions (i.e. CPT and resting state) were investigated using Wilcoxon non-parametric tests for paired samples. The use of nonparametric tests was justified by the non-gaussian distribution of samples, as MIC may assume only positive values in the [0 -1] range.

The threshold for statistical significance was corrected in accordance with the Bonferroni rule considering a total of 10 multiple comparisons (5 EEG bands × 2 MF cumulants), thus α = α 1 /10 = 0.005, with α = 0.05 the uncorrected significance threshold.

III. EXPERIMENTAL RESULTS

Figures 2 and3 show the topographic distribution of MIC differences between CPT and resting phase for c 1 and c 2 , respectively. Particularly, while the c 1 parameter show nonspecific patterns among channels with positive and negative values in the [-0.5, 0.5] range, c 2 shows MIC values within a [-0.2, 0.2] range with a prevalence for negative values (i.e., higher MIC values are associated with resting state), especially in the β and γ bands. Consequently, the functional relationship between brain and heart MF cumulants is higher during resting state than CPT.

Figure 4 shows the outcome of statistical comparisons between MIC values derived from resting state vs. CPT expressed in terms of -log 10 (P -value). The significance threshold is set to -log 10 (α) = -log 10 (0.005) = 2.301 because of the Bonferroni correction. While no statistically significant comparisons are associated with the first-order MF cumulant c 1 , several significant differences were found at all EEG frequency bands for c 2 .

IV. DISCUSSION AND CONCLUSION

We reported on a novel methodology for quantifying brainheart interplay based on functional relationships between MF cumulants. The method was tested on experimental EEG and heartbeat series gathered from 24 healthy volunteers undergoing a CPT elicitation that follows a resting state session. Linear and nonlinear features were derived from the cumulants of the MF spectrum of brain and heartbeat series similarly to [START_REF] Catrambone | Heartbeat dynamics analysis under coldpressure test using wavelet p-leader non-gaussian multiscale expansions[END_REF], [START_REF] Catrambone | Wavelet p-leader non-gaussian multiscale expansions for eeg series: an exploratory study on cold-pressor test[END_REF], and their concurrent variation across time-scales was estimated through the MIC. Results suggest that significant statistical differences between the two experimental conditions may be associated with the nonlinear MF feature c 2 , and no differences may be linked to the linear quantifier c 1 . Higher c 2 values are in the resting state session with respect to the CPT session.

Previous studies investigated MF properties of EEG and heartbeat series separately during CPT and showed that the stressful condition leads to a decrease in MF properties towards a more linear dynamics [START_REF] Catrambone | Heartbeat dynamics analysis under coldpressure test using wavelet p-leader non-gaussian multiscale expansions[END_REF], [START_REF] Catrambone | Wavelet p-leader non-gaussian multiscale expansions for eeg series: an exploratory study on cold-pressor test[END_REF]. Here, we found that the functional coupling between brain-and heartbeatderived MF metrics, quantified through MIC, follow the same trend as the single systems. These results confirm lower brain's and heart's multifractality during CPT than resting state, and suggest that a consistent decrease in their functional coupling exists in case of a CPT elicitation. On the other hand, linear properties of the MF spectrum (c 1 ) seem to be associated with a high inter-subject variability, therefore minimizing eventual statistical difference between experimental sessions. Note that p-values associated with the c 2 cumulant are as low as 10 -8 , and more pronounced differences between the two experimental conditions are associated with higher EEG oscillations in the β and γ bands.

The proposed methodology presents a limitation related to the non-quantification of the directionality of the functional coupling between brain and heartbeat dynamics. Moreover, statistical analyses were performed at a channel-wise level, thus limiting any inference regarding specific brain regions.

This study represents a first attempt for the quantification of linear and nonlinear dependence between MF properties of brain and heartbeat dynamics. We conclude that a functional nonlinear relationship between brain-and heartbeat-related MF cumulants exist, and higher values are associated with the resting state than the CPT. Further endeavours will be directed to the investigation of other features derived from the MF spectrum, as well as to the application of physiological data collected during different experimental conditions.

Fig. 1 :

 1 Fig. 1: Schematic representation of the implemented processing pipeline. p-leaders (p)

Fig. 4 :

 4 Fig. 4: Negative logarithm 10 of the p-values from Wilcoxon non-parametric tests for the resting state vs. CPT comparison using MIC values at all EEG channels. Rows indicate MF cumulants, and columns indicate EEG-PSD frequency bands vs. µ RR . Significant values are higher than -log 10 (α = 0.005) = 2.301.

Bioengineering and Robotics Research Center E.Piaggio & Department of Information Engineering, School of Engineering, University of Pisa, Italy;

IRIT-ENSEEIHT, Université de Toulouse, CNRS, France;

Department of Electronics, Informatics and Bioengineering, Politecnico di Milano, Italy;

Univ Lyon, ENS de Lyon, Univ Claude Bernard, CNRS, Laboratoire de Physique, Lyon, France; * Corresponding author. Email at: vincenzo.catrambone@ing.unipi.it The research has received partial funding by the Italian Ministry of Education and Research (MIUR) in the framework of the CrossLab project (Departments of Excellence).


Dr Petina L Pert
email: petina.pert@csiro.ausuggested

Lauric Thiault

Matthew I Curnock

Susanne Becken

Joachim Claudet

Beauty and the reef: Evaluating the use of non-expert ratings for monitoring aesthetic values of coral reefs

Keywords: Aesethetic, value, beauty, appreciation, cultural ecosystem services, attitudes, perception, Great Barrier Reef

UNESCO World Heritage sites are affected by global change. The Great Barrier Reef (GBR) is an iconic World Heritage site whose values for which the site was inscribed, including its natural beauty, are dramatically declining. However, there is currently no program or documented methods for long term monitoring of aesthetics for coral reefs. Here, Here, we quantified and assessed the drivers of the GBR aesthetic value perception form a survey of 1,417 individuals who rated 181 photos on a 10-point scale for their level of attractiveness. Each photo was rated at least 380 times. Coral topography, fish abundance and visibility were positively correlated 1 with aesthetic ratings. Older people and those with interest in coral reefs and who dived on the reef rated photos more highly. Reporting on aesthetic values of the UNESCO Whorld Heritage Sites will become increasingly as declines in their ecological state continue.

Introduction

The reframing of nature by the Millennium Ecosystem Assessment report (Millennium Ecosystem Assessment, 2003), in particular, as providing goods and services, has enabled a mechanism by which the importance of natural places can be reconsidered [START_REF] Diaz | Linking functional diversity and social actor strategies in a framework for interdisciplinary analysis of nature's benefits to society[END_REF].

The aesthetic services that are provided by nature have become formally acknowledged and valued, where the provision of aesthetic services is strongly related to all other cultural services including recreation [START_REF] Ghermandi | Recreational, Cultural, and Aesthetic Services from Estuarine and Coastal Ecosystems[END_REF][START_REF] Junge | Aesthetic quality of agricultural landscape elements in different seasonal stages in Switzerland[END_REF]. Areas with conditions for recreation and leisure are often valuable due to their aesthetic and cultural heritage and identity qualities (e.g. [START_REF] Chan | Where are Cultural and Social in Ecosystem Services? A Framework for Constructive Engagement[END_REF][START_REF] Comberti | Ecosystem services or services to ecosystems? Valuing cultivation and reciprocal relationships between humans and ecosystems[END_REF][START_REF] Gee | Cultural ecosystem services in the context of offshore wind farming: A case study from the west coast of Schleswig-Holstein[END_REF][START_REF] Jobstvogt | Looking below the surface: The cultural ecosystem service values of UK marine protected areas (MPAs)[END_REF]). Likewise, areas with scenic beauty often provide inspiration and opportunities for education, holding significant existence values (e.g., [START_REF] De Oliveira | What value São Pedro's procession? Ecosystem services from local people's perceptions[END_REF][START_REF] Fletcher | Revealing marine cultural ecosystem services in the Black Sea[END_REF][START_REF] Hashimoto | Mapping and characterizing ecosystem services of social-ecological production landscapes: case study of Noto, Japan[END_REF][START_REF] Martínez-Pastur | Spatial patterns of cultural ecosystem services provision in Southern Patagonia[END_REF][START_REF] Outeiro | Framing local ecological knowledge to value marine ecosystem services for the customary sea tenure of aboriginal communities in southern Chile[END_REF]). Furthermore, De Groot et al. (2002) refers to 'aesthetic information' as an ecosystem service category, provide 'scenery' as an example of an aesthetic service, and a feeling of 'enjoyment' as the benefit of this service. [START_REF] Boyd | What are ecosystem services? The need for standardized environmental accounting units[END_REF] suggest that 'amenity and fulfilment' are examples of aesthetic benefits and 'natural land cover in viewsheds' is the ecosystem service providing the benefits. They argue this differentiation is essential in defining a workable accounting system to facilitate land-use decision-making and management.

For the purpose of this paper, the term aesthetic value(s) has been taken as defined by Context Pty Ltd (2013 p. 33), being: ".. the response (the aesthetic response) derived from the experience of an environment or parts of an environment. Human sensessight, touch, smell, sound, movement -are important in how humans experience an environment. And culture, knowledge, expectations and experience mediate sensory perceptions. Aesthetic response can therefore be said to be linked to:

 the characteristics of an environment  culturally or personally derived preferences." Furthermore, we have summarised the by the many contributing environmental and experiential characteristics that contribute to aesthetic response, satisfaction and derived value in a natural setting (Fig. 1). The natural beauty and aesthetic appreciation of landscapes has previously been studied in various terrestrial ecosystems around the world [START_REF] Arriaza | Assessing the visual quality of rural landscapes[END_REF][START_REF] Beza | The aesthetic value of a mountain landscape: a study of the Mt. Everest Trek[END_REF][START_REF] De Pinho | Influence of aesthetic appreciation of wildlife species on attitudes towards their conservation in Kenyan agropastoralist communities[END_REF][START_REF] Frank | Assessment of landscape aesthetics-Validation of a landscape metrics-based assessment by visual estimation of the scenic beauty[END_REF][START_REF] Gobster | The shared landscape: what does aesthetics have to do with ecology?[END_REF][START_REF] Hoffman | Silviculture and forest aesthetics within stands, The New York center for forestry research and development[END_REF][START_REF] Howley | Landscape aesthetics: Assessing the general publics' preferences towards rural landscapes[END_REF][START_REF] Rosley | Perceiving the Aesthetic Value of the Rural Landscape Through Valid Indicators[END_REF][START_REF] Schirpke | Predicting scenic beauty of mountain regions[END_REF][START_REF] Sheppard | Visual analysis of forest landscapes[END_REF][START_REF] Van Den Berg | Group differences in the aesthetic evaluation of nature development plans: a multilevel approach[END_REF][START_REF] Van Der Jagt | Unearthing the picturesque: The validity of the preference matrix as a measure of landscape aesthetics[END_REF][START_REF] Van Zanten | Continental-scale quantification of landscape values using social media data[END_REF]Wet Tropics Management Authority, 2017;[START_REF] Xu | Valuing Biodiversity, Aesethetics, and Job Losses Associated with Ecosystem Management Using Stated Preferences[END_REF], however very few studies have occurred in marine ecosystems [START_REF] Fenton | Perception and evaluation of the coastal zone: implications for coastal zone planning[END_REF] or more specifically within coral reefs [START_REF] Dinsdale | Assessing coral reef condition: eliciting community meanings[END_REF][START_REF] Fenton | Re-presenting the great barrier reef to tourists: Implications for tourist experience and evaluation of coral reef environments[END_REF]. Coral reef ecosystems worldwide are significant sources of consumptive and non-consumptive economic values [START_REF] Wielgus | Dose-Response Modeling of Recreationally Important Coral Reef Attributes: A Review and Potential Application to the Economic Valuation of Damage[END_REF] and supply vast numbers of people with goods and services such as seafood, tourism and recreational opportunities, coastal protection, as well as aesthetic and cultural benefits [START_REF] Martin | A systematic quantitative review of coastal and marine cultural ecosystem services: Current status and future research[END_REF][START_REF] Moberg | Ecological goods and services of coral reef ecosystems[END_REF][START_REF] Werner | Rapid prioritization of marine ecosystem services and ecosystem indicators[END_REF]. Many attributes have been described as to what contributes to these non-consumptive aesthetic values including the abundance of fish and corals [START_REF] Dixon | Meeting Ecological and Economic Goals: Marine Parks in the Caribbean[END_REF][START_REF] Wilkinson | Global Change and Coral Reefs: Impacts on Reefs, Economies and Human Cultures[END_REF][START_REF] Williams | Differences between Protected and Unprotected Reefs of the Western Caribbean in Attributes Preferred by Dive Tourists[END_REF], existence of iconic species/habitats [START_REF] Jobstvogt | Looking below the surface: The cultural ecosystem service values of UK marine protected areas (MPAs)[END_REF][START_REF] Levin | Developing meaningful marine ecosystem indicators in the face of a changing climate[END_REF][START_REF] Werner | Rapid prioritization of marine ecosystem services and ecosystem indicators[END_REF], and water clarity [START_REF] Dixon | Meeting Ecological and Economic Goals: Marine Parks in the Caribbean[END_REF][START_REF] Hatton Macdonald | Valuing coastal water quality: Adelaide, South Australia metropolitan area[END_REF][START_REF] Lee | Appearance's Aesthetic Appreciation to Inform Water Quality Management of Waterscapes[END_REF][START_REF] Tosic | Selection of parameters in the design of beach coastal water quality monitoring programs[END_REF][START_REF] Wilkinson | Global Change and Coral Reefs: Impacts on Reefs, Economies and Human Cultures[END_REF].

The development of indicators and quantifying the aesthetic value that people place on different landscapes, seascapes, and ecosystems has been approached in a number of ways by researchers over time. Primarily these have focused on using landscape assessment methods which either focus on understanding experiential preferences, the physical attributes/attractiveness of a place [START_REF] Tardieu | From supply to demand factors: What are the determinants of attractiveness for outdoor recreation?[END_REF], landscape preferences, [START_REF] Atauri | Visitors' landscape preferences as a tool for management of recreational use in natural areas: a case study in Sierra de Guadarrama (Madrid, Spain)[END_REF][START_REF] De La Fuente De Val | Relationship between landscape visual attributes and spatial pattern indices: a test study in Mediterranean-climate landscapes[END_REF][START_REF] De Lucio | Landscape preferences and behavior of visitors to Spanish National Parks[END_REF][START_REF] Huang | Visitor responses to the changing character of the visual landscape as an agrarian area becomes a tourist destination: Yilan County, Taiwan[END_REF][START_REF] Múgica | The role of on-site experience on landscape preferences. A case study at Doñana National Park (Spain)[END_REF], scenic quality, or the experience of the place. Two different methodological approaches that assess landscape aesthetics either from an objective or a subjective point of view have been established in the past [START_REF] Frank | Assessment of landscape aesthetics-Validation of a landscape metrics-based assessment by visual estimation of the scenic beauty[END_REF]. More recent approaches such as those by [START_REF] Casalegno | Spatial covariance between aesthetic value & other ecosystem services[END_REF] and [START_REF] Van Zanten | Continental-scale quantification of landscape values using social media data[END_REF] have used crowd-source datasets, machine learning, [START_REF] Haas | Can we measure beauty? Computational evaluation of coral reef aesthetics[END_REF][START_REF] Seresinhe | Quantifying the Impact of Scenic Environments on Health[END_REF][START_REF] Spalding | Mapping the global value and distribution of coral reef tourism[END_REF] and geo-tagged digital photographs uploaded to social media platforms to quantify and objectively measure the perceived aesthetic value of ecosystems.

The emergence of aesthetic values methods in the Australian heritage practice although only new, is important as there is no established World Heritage methodology for assessing natural beauty and aesthetic values under Criterion vii (i.e. contain superlative natural phenomena or areas of exceptional natural beauty and aesthetic importance). In the Great Barrier Reef and Wet Tropics World Heritage Areas in Australia attempts have been made to define and measure aesthetic values (Context Pty Ltd, 2013;Wet Tropics Management Authority, 2017). More recently new methods have been developed utilising non-expert visual assessments [START_REF] Marshall | Identifying indicators of aesthetics in the Great Barrier Reef for the purposes of management[END_REF] and innovative technologies using artificial intelligence [START_REF] Becken | Monitoring aesethetic value of the Great Barrier Reef by using innovative technologies and artificial intelligence[END_REF]. Given the importance of the tourism industry in these areas, those features and attributes that currently draw visitors to these areas are important, particularly to tourism operators in the areas. Furthermore, a number of values for which the Great Barrier Reef was originally inscribed for are on the decline after back to back coral bleaching events recently in 2016 and 2017, terrestrial run-off of sediments and nutrients and associated Crown-of-Thorns Starfish (COTS) (Acanthaster planci) outbreaks, acidification associated with climate change, impacts from severe cyclones, and coral diseases [START_REF] Brodie | A critical review of environmental management of the 'not so Great' Barrier Reef[END_REF][START_REF] De'ath | The 27-year decline of coral cover on the Great Barrier Reef and its causes[END_REF][START_REF] Hughes | Securing the future of the Great Barrier Reef[END_REF][START_REF] Hughes | Global warming and recurrent mass bleaching of corals[END_REF], which have all affected the aesthetics and natural beauty of some parts of the Region above and below the water (Great Barrier Reef Marine Park Authority, 2019).

Appealing features of methods which use non-expert rating scores of aesthetic beauty include simplicity and scalability, and as [START_REF] Marshall | Identifying indicators of aesthetics in the Great Barrier Reef for the purposes of management[END_REF] and [START_REF] Haas | Can we measure beauty? Computational evaluation of coral reef aesthetics[END_REF] have shown, relevance/correlation to environmental attributes associated with Reef health. However, before such a method can be implemented and used in a management context, managers need to understand the method's sensitivity, extrapolative power, and any potential biases associated with non-expert visual ratings of scenes. In this study we evaluate (a) the sensitivity of rating scores used in [START_REF] Marshall | Identifying indicators of aesthetics in the Great Barrier Reef for the purposes of management[END_REF] to compare underwater coral reef scenes, (b) the sampling design requirements to establish appropriate confidence levels, and (c) biases associated with demography and knowledge levels/expertise.

Methods

Study site

The Great Barrier Reef (GBR) of Australia is one of the world's largest coral reef systems. It houses globally-significant biodiversity, a rich diversity of underwater habitats, and a scenic appeal that is universally recognised [START_REF] Hughes | Climate change, human impacts, and the resilience of coral reefs[END_REF]. These values make it an international tourism attraction: in 2013, the Great Barrier Reef received an estimated 53.3M days of use, with 98% of all days comprising tourism visitation [START_REF] Marshall | Advances in monitoring the human dimension of natural resource systems: an example from the Great Barrier Reef[END_REF]. Tourism operators and commercial fishers are particularly appreciative of the natural beauty of the region, with the tourism industry alone worth over $5.2 billion each year [START_REF] Marshall | Advances in monitoring the human dimension of natural resource systems: an example from the Great Barrier Reef[END_REF].

It also holds immense value for local communities, where 'beautiful' is one of the first words that come to local residents' minds when thinking of the GBR [START_REF] Marshall | Advances in monitoring the human dimension of natural resource systems: an example from the Great Barrier Reef[END_REF]. Residents and tourists seek opportunities to appreciate the natural aesthetics of the GBR through experiences such as diving and snorkelling, visiting beaches, boating and scenic flights [START_REF] Marshall | Advances in monitoring the human dimension of natural resource systems: an example from the Great Barrier Reef[END_REF].

Data analysis

We developed a linear mixed effect model to generate expectations of how image scores were related to eight covariables expected to influence the aesthetic score. These included coral size, fish size, fish diversity, presence of algae, coral cover, topography, fish density, and visibility [START_REF] Marshall | Identifying indicators of aesthetics in the Great Barrier Reef for the purposes of management[END_REF]. Respondents and images were used as random effects. Variance inflation factor (VIF) scores 3.28, indicating low collinearity among covariates. We then identified XXXX and YYYY, defined as respondents scoring more than two standard deviations higher or lower than expectations from our model, respectively. Our approach thus reveals the respondents that tend to provide more extreme scores given the general characteristics of the images they evaluate. Finally, we assessed the link between socioeconomic characteristics and whether a respondent was XXXX, average, or YYYY using a Fisher's exact test.

In order to help practitioners find a balance between low cost and high accuracy, we examined how error in the image scores varied according to the number of respondents (i.e. sample size). Specifically, we calculated the error (i.e. score's 95% confidence interval) at varying sample sizes (N = 3 to 380 respondents in increments of 2, taken randomly without replacement from the original pool of respondents). This process was repeated 100 times for each of the 181 images in order to calculate the median error each image add at various sample size.

Five indicators were selected that reflected distinct attributes of coral reefs: (i) coral cover; (ii) coral pattern; (iii) coral topography; (iv) fish abundance; and (v) visibility [START_REF] Marshall | Identifying indicators of aesthetics in the Great Barrier Reef for the purposes of management[END_REF].

All analyses were implemented in R statistical analysis software v3.4.0.

Results

The average aesthetic score ranged between X and X. Four images (Fig. 2) indicate the frequency and average aesthetic score.

Discussion

• CHALLENGES Reporting through monitoring, and actively managing, aesthetic quality is potentially a critically useful strategy for natural resource managers that aim to effectively deliver both conservation and social outcomes. The advantages of monitoring aesthetic values lie essentially around the opportunity to better engage with the public, correspondingly achieving conservation [START_REF] Frank | Assessment of landscape aesthetics-Validation of a landscape metrics-based assessment by visual estimation of the scenic beauty[END_REF]. People appreciate and relate to aesthetics more than any other aspect of the natural environment [START_REF] Marshall | Measuring What Matters in the Great Barrier Reef[END_REF][START_REF] Marshall | Advances in monitoring the human dimension of natural resource systems: an example from the Great Barrier Reef[END_REF][START_REF] Marshall | The Dependency of People on the Great Barrier Reef, Australia[END_REF]. Aesthetics is also a unifying value associated with natural places, whilst biodiversity is often a divisive factor [START_REF] Hill | Reconciling Biodiversity Indicators to Guide Understanding and Action[END_REF][START_REF] Lecuyer | Building on common ground to address biodiversity conflicts and foster collaboration in environmental management[END_REF][START_REF] Schmidt | Biodiversity Conservation and Indigenous Land Management in the Era of Self-Determination[END_REF][START_REF] Young | The emergence of biodiversity conflicts from biodiversity impacts: characteristics and management strategies[END_REF]. In the case of iconic places such as the Great Barrier Reef, aesthetic values were a major factor in the establishment of World Heritage Area status, and the reporting of aesthetic quality and its maintenance is likely to become mandatory. For these reasons, efforts are currently underway to explore how aesthetic quality might be recorded and reported on within the Great Barrier Reef (NESP reference).

Perhaps the reason that aesthetic values have not been explicitly used by natural resource managers to further conservation goals and support is because aesthetics are seen as subjective and essentially unmanageable (Context Pty Ltd, 2013;[START_REF] Dramstad | Relationships between visual landscape preferences and map-based indicators of landscape structure[END_REF]. However, preliminary work by the National Environment and Science Program [START_REF] Becken | Monitoring aesethetic value of the Great Barrier Reef by using innovative technologies and artificial intelligence[END_REF][START_REF] Marshall | Identifying indicators of aesthetics in the Great Barrier Reef for the purposes of management[END_REF] has strongly suggested that developing indicators to rate aesthetic quality is possible [START_REF] Marshall | Identifying indicators of aesthetics in the Great Barrier Reef for the purposes of management[END_REF]. In sum, the development of aesthetic indicators for use in specific environments such as marine environments are still very much in their infancy, but the momentum supporting the development of aesthetic metrics, or indicators, is growing [START_REF] Belhassen | Evaluating the attractiveness and effectiveness of artificial coral reefs as a recreational ecosystem service[END_REF][START_REF] Haas | Can we measure beauty? Computational evaluation of coral reef aesthetics[END_REF][START_REF] Hein | Coral restoration: Socio-ecological perspectives of benefits and limitations[END_REF][START_REF] Jaap | Coral reef restoration[END_REF][START_REF] Tamayo | Estimates of Values of Philippine Reefs' Ecosystem Services[END_REF][START_REF] Turek | Science and technology needs for marine fishery habitat restoration[END_REF][START_REF] Vercelloni | Using virtual reality to estimate aesthetic values of coral reefs[END_REF][START_REF] Vlami | Cultural landscapes and attributes of "culturalness" in protected areas: An exploratory assessment in Greece[END_REF][START_REF] Wessel | An evaluative tool for rapid assessment of derelict vessel effects on coastal resources[END_REF].

The emerging social-ecological field continues to expand methodologies to improve techniques in monitoring aesthetic heritage values, using potential indicators and computations of aesthetic value. However, ongoing examination of which locations or biophysical elements are the most important to the Reef's spectacular seascapes and scenery, remains an information gap.

Conclusion

There are numerous potential management (and other) uses of aesthetic assessments and monitoring with information needs that vary substantially. The rapid assessment approach is potentially useful for some purpose with statistical design requirements and limitations now better understood. Our findings can help World Heritage site managers better contextualise the results and limitations from such monitoring and ensure that future aesthetic monitoring programs can be designed to meet management and reporting needs.

Funding

This research was funded by the National Environmental Science Program (Australia).

Fig. 1 .

 1 Fig. 1. Characteristics and attributes that contribute to aesthetic response, satisfaction and

Fig

  Fig. 2 XXXXXXXXXXXXXXXXXXX

Fig. 4 .

 4 Fig. 4. Sample size describes how many people you would need to assess a photo/coral reef


Benjamin Ribba

Olivier Saut

Thierry Colin

Didier Bresch

Emmanuel Grenier

Jean-Pierre Boissel

A multiscale mathematical model of avascular tumor growth to investigate the therapeutic benefit of anti-invasive agents

Keywords: Inhibitors of Metalloproteinases, Cytostatic agents, Mathematical modeling, Therapeutic development

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

A multiscale mathematical model of avascular tumor growth to investigate the therapeutic benefit of anti-invasive agent

Mathematical models of cancer have been extensively developed with the aim to predict tumor growth and therapeutic strategies efficacy (1; 2; 3; 4; 5). More complex multiscale models have been set up to address questions in the prediction of efficacy and toxicity of chemotherapy, radiotherapy (6; 7), and to investigate the effect of cytostatic treatments such as anti-angiogenic and anti-maturation therapies [START_REF] Arakelyan | A computer algorithm describing the process of vessel formation and maturation, and its use for predicting the effects of anti-angiogenic and anti-maturation therapy on vascular tumor growth[END_REF].

In this paper, we're interested in cytostatic agents called the inhibitors of metalloproteinases (MMPi), developed to slow down cancer growth and to reduce the risk of metastasis. These agents have been designed to inhibit metalloproteinases (MMP), a group of enzymes which are known to play a major role in the degradation of basal membranes and extra-cellular matrixes and thus in cancer invasion. Despite the promising results obtained in animals models, the clinical development of these agents was rather disappointing.

Retrospective analysis have highlighted a couple of potential reasons to explain this failure. Nonetheless, other studies are required to allow the resumption of further clinical trials.

The aim of this paper is to present a mathematical model which integrates basic understanding on the role of MMP in cancer growth in order to provide some explanation on the therapeutic failure and prospectively to help the design of further clinical trials of MMPi. Our model is composed by an agestructured cell cycle model with three populations. We applied Darcy's law to simulate spatially tumor dynamic. The efficacy of simulated MMPi treatment has been evaluated in term of tumor growth reduction. [START_REF] Kufe | Cancer Medicine[END_REF]. Many evidences in experimental studies prior to the 1990 reveal the role of MMP in promoting cancer growth, tissue invasion, and development of metastasis [START_REF] Giavazzi | Preclinical development of metalloproteasis inhibitors in cancer therapy[END_REF]. Despite the structural differences between the members of the MMP family, non-specific inhibitors have been developed since at the time the precise role of each members was unknown and MMP action was seen globally to promote invasion and thus cancer growth [START_REF] Zucker | Critical appraisal of the use of matrix metalloproteinase inhibitors in cancer treatment[END_REF].

Preclinical studies have shown a significant benefit of MMPi to inhibit tumor growth and to prevent metastasis in animal models [START_REF] Coussens | Matrix metalloproteinase inhibitors and cancer: trials and tribulations[END_REF]. Despite these encouraging results, clinical trials have been disappointing. Phase I showed toxicity which was not expected from preclinical studies. Results of Phase II were difficult to analyze [START_REF] Coussens | Matrix metalloproteinase inhibitors and cancer: trials and tribulations[END_REF]. Indeed, the conventional endpoints used to measure efficacy of anti-cancer drugs, e.g., reduction of tumor volume, were not relevant for cytostatic agents since they do not kill cell but rather block proliferation. Phase III trials did not show any significant efficacy on the clinical endpoints (14; 15).

Retrospective analysis provided few arguments to explain the failure of these therapeutic developments. It is admitted that a better understanding of each

We present a multiscale mathematical model of avascular tumor growth to investigate the role of MMPi treatment on cancer growth. The mathematical core is composed by an age-structured model where we distinguished between proliferative and quiescent cells. The cell cycle model is embedded in a macroscopic model based on Darcy's law to describe tumor spatial dynamic.

Cell cycle model

We developped an age-structured mathematical model of the cell cycle where cell cycle phase duration were set according to literature [START_REF] Potten | Measure-ment of in vivo proliferation in human colorectal mucosa using bromodeoxyuridine[END_REF]. We distinguished between the proliferative phases P 1 and P 2 , and the quiescent phase Q. We modeled the 'Restriction point' R (17) at the end of phase P 1 where environmental conditions are checked [START_REF] Kufe | Cancer Medicine[END_REF].

Cells evolve in time t and in age a in the proliferative phases of the cell cycle. We assumed that if the local number of cancer cells was above a threshold T h t , then overpopulation was declared and leads cells at the end of P 1 to compartment Q (quiescence) which is not age-structured. Quiescent cells can come back to the proliferative cycle (in phase P 2 ) as soon as the environmental conditions become appropriate again, i.e. the local number of cells goes below T h t . The age-structured model can be written as follows:

             ∂P 1 ∂t + ∂P 1 ∂a = 0, ∂P 2 ∂t + ∂P 2 ∂a = 0, ∂Q dt = (1 -f (N ))P 1 (a = a max,P 1 ) -[f (N )] + Q, (1)
with the following initial conditions:

     P 1 (a = 0) = 2P 2 (a = a max,P 2 ), P 2 (a = 0) = f (N )P 1 (a = a max,P 1 ) + [f (N )] + Q, ( 2
)
where N denotes the local number of cells and f (N ) is a boolean function such that:

f (N ) =      1 if N < T h t 0 otherwise (3)
See Figure 1 for a schematic representation of the cell cycle model.

Macroscopic model

We work in the framework of continuous meachnaics and we descrive the medium by using a velocity v and a pressure p that are related to the Darcy's law (18; 19; 20):

v = -k∇p (4)
where k denotes the media permeability. We refer to [START_REF] Ambrosi | On the closure of mass balance models for tumor growth[END_REF] for a precise description and justification of such models in the context of potential flows. Here, both quantities v and p are related to a mixture involving several species. These species are described by n ϕ (x, y, a, t) that represents the cell numbers (by volume unit) of cancer cells with spatial position (x, y) at time t with age a in a given proliferative cycle phase ϕ. Applying the principle of mass balance, the advection equations are:

     ∂nϕ ∂t + ∂nϕ ∂a + ∇ • (vn ϕ ) = 0, ∀ϕ ∈ {P 1 , P 2 }, dn Q dt + ∇ • (vn Q ) = (1 -f (N ))P 1 (a = a max,P 1 ) -[f (N )] + n Q , ( 5
)
where N (x, y, t) = n Q (x, y, t) + ϕ a n ϕ (x, y) ∀ϕ ∈ {P 1 , P 2 } is the total number of cancer cells (by volume unit). Note that the quiescent phase Q is not age-structured.

We denote M the density of healthy tissue which evolves as follows:

dM dt + ∇ • (vM ) = 0, (6)
and assume that N (x, y, t) + M (x, y, t) = Constant ∀(x, y, t).

Integrating on ages a equations 5, adding them and using the relation N + M is constant, one can find the following relation:

[n ϕ ] a=a max,ϕ a=0 + ∇ • v = (1 -f (N ))P 1 (a = a max,P 1 ) -[f (N )] + n Q . ( 7
)
Now using the initial conditions 2, we find the following equation for the velocity:

∇ • v = P 2 (a = a max,P 2 ). (8)
Thus, the pressure field p satisfies the elliptic equation:

-∇ • (k∇p) = P 2 (a = a max,P 2 ). (9)
We set n ϕ and n Q on part of the boundary where v • ν < 0, with ν denoting the outgoing normal to the boundary . For the pressure p, we set p = 0 on the boundary of the computational domain.

MMP, membrane degradation and action of therapeutic

It has been basically assumed that MMP are secreted by cancer cells and that the degradation of the matrix components was function of MMP. We modeled this as follows:

     M M P = κ M M P N (t), dECM dt = -β(M M P (t)), (10)
where ECM stands for the basal membrane to be degraded by cancer cells and β is an hyperbole-like function written as an 'Emax model':

β(M M P (t)) = β max • M M P (t) κ β + M M P (t) . ( 11
)
where β max and k β are positive constants.

To simulate the evolution of a population of cancer cells under MMPi treatment we basically set M M P (t) = 0 during the whole simulation which directly results in preventing ECM degradation and thus cancer cells dissemination in the surrounding tissue.

Evaluation of the therapeutic benefit

As the action of MMPi is aimed at promoting the passage of cells from proliferation to quiescence, we'll use the proportion of quiescent cells as the efficacy marker. In the following, n Q and n Q will denote respectively the local number of quiescent cells without MMPi and with MMPi treatment. The absolute effect (AE) of the therapy will be compute as follows:

AE = 1 T t=T t=0 ( n Q (t)dxdy N (t)dxdy - n Q (t)dxdy N (t)dxdy ), ( 12
)
where N is the total number of cells when MMPi treatment is applied, and T the global simulation time. The evaluation of the therapeutic benefit will be done as follows:

    
if AE > 0 the treatment is effective otherwise the treatment is useless (or harmful) (13)

Meaningful model parameters

We identified two meaningful parameters which characterize both cell and tissue properties. At the cell level, the threshold T h t is used to decide wether cells continue to proliferate or go to quiescence. Thus, an increase in this threshold characterizes a lack of sensitivity towards overpopulation. At the macroscopic level, the critical parameter is the tissue porosity noted k in the Darcy's law (see Eq. 4). Different values of k results in different cell spatial motion and in consequence different distribution in the cell cycle phases according to the regulation of the cell cycle (see Eq. 3).

Computation domain and initial conditions

We assumed that the domain is composed by healthy and cancer cells delimited from surrounding tissue by a membrane. Figures 2,3, 4 present some illustrations of spatial-temporal evolution of cancer cells in lifelike geometries.

In the following we will keep the spatial configuration presented Figure 4 and note k 0 and k 1 the porosity constants of the healthy tissue and the membrane respectively with k 0 > k 1 (see Figure 5). Initially, cells are distributed uniformly in the phases of the cell cycle.

Parameters used in the simulations are presented Table 1.

Results

According to the macroscopic model (Darcy's law), the media porosity regulates the tendency for cells to move towards new areas or to remain at their spatial position. When the MMPi treatment is applied, it prevents membrane degradation. The porosity of the membrane should not allow cells to invade new tissue areas. In consequence, proliferating cells should remain at their spatial location. After a certain number of cycles, we expect that the local number of cells goes over the overpopulation threshold T h t which promote proliferative cells to move to quiescence. This will result in a positive absolute effect (see Eq. 12).

Influence of porosity ratio k 0 /k 1 on the efficacy of MMPi treatment

We first analyze how the efficacy of MMPi treatment is sensitive to the variation of the porosity ratio k 0 /k 1 . We ran six simulations for ratio porosity values 1, 10 2 , 10 4 , 10 8 , 10 10 , 10 12 . For each porosity ratio values, we ran the simulations over 300 time steps, with and without MMPi treatment, and compute the absolute benefit with Eq. 12. Results are presented Figure 6.

While the ratio k 0 /k 1 increases (the porosity of the membrane decreases respect the porosity of the surrounding tissue), we get as expected a positive absolute effect which indicates that the MMPi treatment is effective, i.e., it promotes proliferative cells to move massively to quiescence. Note that when the porosity of the surrounding tissue is equal to the membrane porosity, there's no therapeutic benefit to apply the MMPi treatment. Indeed, when the tissue porosity is equal to the membrane's one (k 0 /k 1 = 1), degradation does not have any effect on spatial availability for cells and thus repartition in phase P and Q of the cell cycle is not affected.

Influence of overpopulation threshold (T h t ) on the efficacy of MMPi treatment

We ran six simulations for overpopulation threshold values 1000, 2000, 3000, 4000, 6000, and 8000. For each of these values, we ran the simulations over 300 time steps, with and without MMPi treatment, and compute the absolute benefit with Eq. 12. For relatively low value of T h t we get a positive absolute effect which means that MMPi treatment is effective in slowing down tumor growth. When T h t is increased beyond a critical value (around 4000), the absolute effect collapses and becomes negative. Results are displayed Figure 7 Variation of the overpopulation threshold can change drastically the effect of the MMPi treatment: from beneficial to harmful.

We plotted Figure 8 the proportion of quiescent cells over time with and without MMPi for three different overpopulaton threshold values (1000,4000, and 8000). For the two first values of the threshold (1000, and 4000) the absolute benefit is positive (see Figure 7) while it is negative for the ultimate value (8000). Indeed, for T h t = 8000, quiescent cells in the case of MMPi treatment appear with a delay respect the situation without MMPi. The negative absolute effect results essentially from this delay. We show in Figure 9 the spatial configuration of quiescent and proliferative cells without MMPi (top) and with MMPi (bottom) from three different time (from left to right). The difference in the evolution of cell spatial configuration with and without MMPi can explain the delay observed in the appearance of quiescent cells when no MMPi is applied. When simulated the tumor growth with the presence of the treatment (bottom), tumor cells tend to extend uniformly within the computational domain which generate a homogeneous spheroid-like shape. On the contrary, in the absence of treatment, membrane degradation occurs and cells have tendency to be engulfed in the area when the membrane has been degraded. Cells at the center and the top of the initial tumor mass will be attracted towards the hole in the membrane and will not expand uniformly has before. In this case, the central part of the tumor mass will become quiescent more quickly than in the case when membrane degradation is prevented.

Conclusions

MMPi were thought as promising therapeutic agents to slow down cancer growth and reduce the risk of metastasis. Their clinical development have failed despite good results obtained in preclinical trials. Retrospective analysis have highlighted the lack of knowledge on the role of MMP in cancer progression and several studies have emphasized the need to apply MMPi treatments on early stage cancer patients. It is known nowadays that MMP can also, according to the family enzymes they belong to and/or according to the patient disease stage, inhibit cancer growth [START_REF] Pavlaki | Matrix metalloproteinase inhibitors (MMPIs): the beginning of phase I or the termination of phase III clinical trials[END_REF]. In this paper, we somehow investigated wether these conclusions could have been anticipated by integrating within a multiscale mathematical model the knowledge that was available at the time.

We built a mathematical model which integrates the basic actions of MMP within a macroscopic tumor growth regulated by a continuous cell cycle model. Simulation results obtained show that the overall efficacy of MMPi treatment may depend on the multiscale coupling between cell cycle regulation and tumor macroscopic behavior. Indeed, with the simple spatial configuration used, we showed that as cancer cells become insensitive to anti-growth signal, MMPi treatment efficacy decreases.

However, as shown, the results obtained are sensitive to the spatial configuration of tumor cells and other initial configuration taken may lead to different and perhaps opposite results.

Nevertheless, the presented results are comparable to the recent conclusions raised by the retrospective analysis. Indeed advanced cancer are usually characterized by several mutations which give to cells a proliferating advantage by loosing their sensitivity to anti-growth signals. The simulation results emphasize the role of cells loosing sensitivity to anti-growth signals in the lack of efficacy of MMPi treatment. This may explain the lack of efficacy obtained during clinical development where MMPi treatment was given only at advanced stage cancer patients. However, the MMPi action preventing membrane degradation and tumor cells to invade the surrounding tissue should have been taken into account, together with the effect on cell proliferation, for a proper evaluation of therapy efficacy.

On a technical point of view, a three-dimensional tumor growth model could reveal new factors in the dynamics. We chose a continuous approach which provides density cells rather than actual cell number and we didn't consider cell shape, which has been shown to be important in the description of growth control process [START_REF] Folkman | Role of cell shape in growth control[END_REF]. To be more realistic, the present model in its 3D version should be coupled with an angiogenic model to allow simulation of MMPi treatment action on vascular tumor growth.

In conclusion, these results should be taken as theoretical and qualitative. Nonetheless, they highlight the need to consider multiscale regulation of tumor growth in the prediction of cancer treatment efficacy. Beyond the complexity of the molecular machinery, prediction of cytostatic therapies efficacy must integrate multiscale regulation of tumor growth via the cell cycle dynamic.

Nowadays, efforts are made to develop specific MMPi inhibitors to target the MMP which are effectively responsible for membrane degradation. Hopefully, a mathematical model as the one presented here may be used to optimize the design of further clinical development of these agents.

  Ribba, Olivier Saut, Didier Bresch, Emmanuel Grenier, Jean-Pierre Boissel To cite this version: Benjamin Ribba, Olivier Saut, Didier Bresch, Emmanuel Grenier, Jean-Pierre Boissel. A multiscale mathematical model of avascular tumor growth to investigate the therapeutic benefit of anti-invasive agent. Journal of Theoretical Biology, 2006, 243 (4), pp.532-541. 10.1016/j.jtbi.2006.07.013. hal-1 Introduction

Fig. 1 .

 1 Fig. 1. The proliferative cycle is composed by two phases: P 1 and P 2 . At the end of P 2 , cells divide and new born cells begin their cycle in P 1 . At the last age of phase P 1 , we modeled the restriction point R, where external conditions are checked. If overpopulation occurs, proliferative cells go to the Q compartment (quiescence). These cells can go back in the proliferative cycle (at the first step of phase P 2 ) if the external conditions allow it.

Fig. 2 .

 2 Fig. 2. From left to right and top to down. Illustration of ECM degradation and tumor cells in vitro invasion. Initially, tumor cells compose a central spot (red then light) surrounded by ECM (dark red). The MMP degrade the ECM (black ring) and allow tumor cells to invade the domain.

Fig. 3 .

 3 Fig. 3. From left to right and top to down. Illustration of in situ ductal carcinoma growth without basal membrane degradation. Initially, tumor cells compose a small area located at the bottom-left area of the epithelial layer. Tumor cells proliferate and deform the basal membrane (white ring) of the duct.

Fig. 4 .

 4 Fig. 4. From left to right and top to down. Illustration of invading epithelial cancer cells with basal membrane and ECM degradation. Initially, tumor cells compose a central small area adjacent to the membrane. Tumor cells proliferate, MMP degrade the membrane and cells invade the surrounding tissue.

Fig. 5 .

 5 Fig.5. Configuration used for the computer simulations aimed at evaluating the therapeutic benefit of MMPi treatment. Cancer cells, which form initially a small square, are embedded in the healthy tissue with porosity k 0 and placed against a basal membrane with higher porosity k 1 .

Fig. 6 .

 6 Fig. 6. Sensitivity of the absolute effect (AE) to the variation of the porosity ratio between surrounding tissue and basal membrane (k 0 /k 1 ).

Fig. 7 .

 7 Fig. 7. Sensitivity of the absolute effect (AE) to variation of the overpopulation signal (T h t ).

Fig. 8 .

 8 Fig. 8. Proportion of quiescent cells (quiescent rate) over simulation time (in hour) without MMPi (black curves) and with MMPi (red curves) for three overpopulation threshold values: T h t = 1000, 4000, 8000.

Fig. 9 .

 9 Fig. 9. From left to right: Spatial configuration of proliferative cells (white), and quiescent cells (red, composing the inner part of the mass) at time t = 40, 128, and 300. Top: Without MMPi. Bottom: With MMPi treatment. At time t = 128, quiescent cells have appeared in the case of no MMPi treatment while at the same time, all cells are proliferative in the case of MMPi treatment.

  In vivo tumor growth is often predicted by the Gompertz model (1; 2) which features that at the beginning of the growth process, a large part of the cells are cycling in the proliferative phases which results in tumor exponential growth. As the number of cells increases, the local environmental can not supply cells with enough space, nutrient, and oxygen, and some cells are moving to quiescence. During cancer growth, new genetic mutations confer to the cells the ability to become unsensitive to anti-growth signals and thus to continue to proliferate.

	Anti-neoplastic agents can be divided in two categories. The first, called cyto-
	toxics, are aimed at killing cells by damaging their genetic contents. Standard
	chemotherapies and ionizing radiations belong to this class. Often their use is
	limited by the toxicity they generate on healthy tissue. Recently, the better
	understanding of the molecular circuits of cancer (10) encouraged the develop-
	ment of targeting treatments acting on molecular processes involved in tumor
	cancer in order to block cancer growth without damaging healthy cells. In
	interfering with the molecular pathways involved in cell proliferation (see (9)
	for a summary of these molecular pathways), these cytostatic agents have the
	ability to promote proliferative cells to move to quiescence and dormant cells
	to remain blocked in the quiescent phase of the cell cycle. Examples of such
	therapeutic agents are Cetuximab which prevent Epidermal Growth Factor
	(EGF) to bind EGFR and activate the mitogenic cascade (MAPK). The anti-
	angiogenic agents Bevacizumab prevent the binding of the Vascular Endothe-
	lial Growth Factor (VEGF) to the receptors FLT-1 and FLK-1 and inhibit
	angiogenesis. In this paper, we are interested in anti-invasive cytostatic agents
	called the inhibitors of metalloproteinases (MMPi). Matrix Metalloproteinases
	(MMP) are a group of enzymes able to degrade the components of the Extra
	Cellular Matrix (ECM)
	During normal morphogenesis, cells evolve in a cell cycle which we often de-
	composed by a proliferative cycle where cells evolve to duplicate DNA and
	divide, and a quiescent phase where cells rest and wait for proliferative signals.
	It is known that proliferative cells can go to quiescence in case of anti-growth
	signals such as hypoxia, i.e. lack of oxygen, and/or nutrient deprivation. How-

2 Biological background

Solid cancers are characterized by an excessive proliferation of cells from epithelial layers. Due to several genetic mutations, cancer cells acquire the ability to produce growth signals and to loose sensitivity to anti-growth signals; the ability to escape from death processes, to replicate indefinitely, and to invade surrounding tissue to form metastasis (see

[START_REF] Hanahan | The hallmarks of cancer[END_REF]

for further details). ever, quiescent cells can turn back in the proliferative cycle if growth signals arise.

  .

	Parameter	Description	Equation	Value	Source
	T	Global simulation time
	a max,P 1	Duration of P 1	1 and 2	20	(16)
	a max,P 2	Duration of P 2	2	18	(16)
	κ M M P	MMP secretion coefficient	10	5 • 10 -4
	κ β	parameter of function β	11	0.01
	β max	parameter of function β	11	0.1
	Table 1
	Dimensionless parameter values used for the numerical simulations.

Acknowledgements

BR is founded by the ETOILE project: 'Espace de Traitement Oncologique par Ions Légers dans le cadre Européen'. We thank Mariama Diallo for reviewing litterature.