Studies concerning the contact heating of food products remain relatively rare in the literature despite the importance of this mode of heat transfer in many industrial and domestic operations such as grilling or pan-frying. This lack of knowledge makes it very difficult to determine rational strategies in order to improve the sensory and nutritional qualities of food products heated by contact.

To deal with this situation, experimental data (water loss, temperature at different positions inside the product) were recorded during contact heating of potato slices performed on a laboratory heating device specially designed for this study. Based on the interpretation of these data, a 2D mathematical model was developed and used to better understand heat and mass transfer phenomena occuring at the contact interface during heating.

VINE

With oil

After 6 seconds of heating (formation of the boiling front)

Contact heating device:

Six operating conditions were tested in order to study the influence of: • The heating power, P heat (3.2 kW.m -2 , 6.4 kW.m -2 and 9.6 kW.m -2 ).

• The presence or not of an oil layer below the product.

Conclusion Experiments

• In the moist region, heat and mass transfer are 2D axisymetric. Heat transfer is due to thermal conduction and water transfer is due to water diffusion. • In the crust, heat and mass transfer are supposed to be z-directional. There is no resistance to vapour transfer and heat transfer is considered to occur in quasisteady state.

Temperature field calculated by the model (COMSOL 5.2):

After 900 seconds of heating (moving of the boiling front)

Thermal contact resistance vs thermal crust resistance:

With oil Without oil

The model was validated by showing good agreements between measured and calculated values of product temperature and water loss. Simulation results show that the overall heating of the product is limited by (i) the boiling front (ii) the developement of the insulating crust. Experimental characterization of the structure of the crust (porosity, permeability) is planned to improve model accuracy.

  Thibaut, FLICK Denis, BERNUAU Emmanuel, BROYART Bertrand UMR Ingénierie Procédés Aliments, AgroParisTech, INRA, Université Paris-Saclay, 91300, Massy, France

Introduction:

Heating protocol: • Pre-heating until the aluminium disc temperature reaches 180°C.

• Product laying and heating at constant heating power, P heat . -𝜆 $? 𝛻𝑇. 𝑛 = ℎ ;E 𝑇 FG (𝑡) -𝑇 Given the high temperature of the aluminium disc, the formation and moving of a boiling front at the bottom surface of the product is assumed. This boiling front divides the product in two regions: the moist region and the dry region (crust). Other assumptions are:

• With oil, the thermal contact resistance (𝑅 EF,;E ) is negligeable compared to the one induced by the crust (𝑅 EF,;< ). Contact heat transfer influences the overall heating of the product only at the early beginning of heating. • Without oil, the thermal contact resistance (𝑅 EF,;E ) is of the same order of magnitude as the thermal resistance of the crust (𝑅 EF,;< ) .


INRODUCTION:

Studies concerning the contact heating of food products remain relatively rare in the literature despite the importance of this mode of heat transfer in many industrial and domestic operations such as grilling or pan-frying. In this work, a 2D mathematical model was developed and used to better understand heat and mass transfer during contact heating of a particular product, the potato.

GOVERNING EQUATIONS:

Given the high temperature of the heating surface, the formation and moving of a boiling front at the bottom surface of the product is assumed. This boiling front is described using a moving interface that divides the product in two regions: the moist region and the dry region (the crust).

RESULTS:

Simulated results are compared to experimental data obtained using the experimental device developed for this study. Results are shown for a constant heating power of 3.2 kW.m -2 6.4 kW.m -2 and 9.6 kW.m -2 with and without the presence of an oil layer below the product.

CONCLUSIONS:

The model was validated by showing good agreements between measured and calculated values of product temperature and water loss. Experimental characterization of the structure of the crust (porosity, permeability) is planned to improve model accuracy (especially for high heating power).

  Aliments, AgroParisTech, INRA, Université Paris-Saclay, 91300, Massy, France

Figure 1 .Figure 3 .

 13 Figure 1. Experimental setup

Figure 4 .

 4 Figure 4. Example of simulation results. Temperature and water content field obtained after 15 min of heating (3.2 kW.m -2 with oil).


Benjamin Ducharne
email: benjamin.ducharne@insa-lyon.fr

Grzegorz Litak
email: g.litak@pollub.pl

Bin Zhang
email: bin.zhang@sdu.edu.cn

Bhaawan Gupta
email: bhaawan.gupta@insa-lyon.fr

Multiscale Characterization and Model for the Dynamic Behavior of Ferroelectric Materials Using Fractional Operators

. Except in memory applications, which are based on polarization switching and hysteresis polarization-electric field relationships, hysteresis is undesired in high-precision sensor, actuator, and capacitor applications. Origins and mechanisms of the piezoelectric hysteresis are complex, and they manifest themselves in qualitatively different forms. An ideal hysteresis loop is

symmetrical, so the absolute value of the positive and negative coercive fields and positive and negative remnant polarizations is equal [START_REF] Boser | Statistical theory of hysteresis in ferroelectric materials[END_REF][START_REF] Morita | Hysteresis model using distribution susceptibility[END_REF]. Hysteresis loops are highly frequency dependent. The main consequences of a frequency increase on hysteresis loops are:

-Increasing coercive fields.

-Decreasing remnant polarizations.

-Varying hysteresis loop areas [START_REF] Potter | Monte carlo simulation of ferroelectric phase transition[END_REF][START_REF] Bartic | Preisach model for the simulation of ferroelectric capacitors[END_REF][START_REF] Sivasubramanian | Equivalent circuit and simulations for the Landau-Khalatniko model of ferroelectric hysteresis[END_REF].

Fractional operators are particularly well adapted to model the frequency dependence of the dielectric hysteresis of a ferroelectric material. Where usual integer derivative operators are usually limited to a relatively weak frequency bandwidth, an approach based on fractional derivatives provides good accuracy between measured hysteresis loops and simulated ones far beyond frequency bandwidth of classical piezoelectric systems. The behavior of piezo ceramics under weak electric fields excitation is usually described by constitutive relations linearized around an operating point. The ceramic comportment as the frequency is varying known as the dielectric relaxation gives important information about the piezo ceramic constitution and about the physical relations describing the polarization behavior. These characteristics are usually monitored using impedance analyzer. Dielectric relaxation is defined as the momentary delay in the dielectric constant of a material. This relaxation is usually described in terms of permittivity as a function of frequency, which can, for ideal systems, be described by the Debye complex formula [START_REF] Debye | Zur Theorie der anomalen dispersion im Gebiete der langwelligen Elektrischen Strahlung[END_REF]:

ε * (ω) = ε ∞ + Δε 1 + iωτ . (10.1)
Here ε ∞ is the sample permittivity under high-frequency excitation, Δε = ε s -ε ∞ where ε s is the quasi-static, low-frequency permittivity, and τ is the characteristic relaxation time. In the case of BaTiO 3 classic piezo ceramic, the Debye equation for the frequency dependence of the complex permittivity is extended to Cole-Cole model [START_REF] Cole | Dispersion and absorption in dielectrics -I alternating current characteristics[END_REF][START_REF] Cole | Dispersion and absorption in dielectrics -II direct current characteristics[END_REF]:

ε * (ω) = ε ′ (ω) -iε ′′ (ω) = ε ∞ + Δε 1 + (iωτ ) α , (10.2)
where ε ′ (ω) and ε ′′ (ω) are the real and imaginary part of the permittivity, respectively. α is linked to the distribution of relaxation time (0 <α<1). Cole-Cole model gives a fractional dependence of the permittivity versus the frequency and gives correct simulation results of a piezo ceramic sample. In case of ferroelectric polymers, an extended version is proposed to consider the unsymmetrical and broadness of the dielectric dispersion curve. This extension called Havriliak-Negami [START_REF] Havriliak | A complex plane representation of dielectric and mechanical relaxation processes in some polymers[END_REF] relaxation requires two fractional parameters:

ε * (ω) = ε ′ (ω) -iε ′′ (ω) = ε ∞ + Δε (1 + (iωτ ) α ) β , (10.3)
with an additional fractional exponent β. In this article, we establish the link between the fractional operators used in the high-amplitude dielectric hysteresis model and the fractional dielectric permittivity models (Cole-Cole), well adapted to the large bandwidth frequency dependence of the permittivity. We also verify if the dynamic mechanisms responsible to the hysteresis frequency dependence are similar to those linked to the frequency dependence of the permittivity.

Nonlinear Model Using a Fractional Derivative

High Excitation Amplitude Fractional Hysteresis Model

The high-amplitude fractional dielectric hysteresis model is used in this study. This model has already been tested successfully previously [START_REF] Ducharne | Low frequency modelling of hysteresis behaviour and dielectric permittivity in ferroelectric ceramics under electric field[END_REF][START_REF] Guyomar | Time fractional derivatives for voltage creep in ferroelectric materials: theory and experiment[END_REF]. Furthermore, it is described precisely with a large number of comparisons of measure/simulation to illustrate its efficiency [START_REF] Zhang | Energy harvesting based on piezoelectric Ericsson cycles in a piezo ceramic material[END_REF][START_REF] Guyomar | High frequency bandwidth polarization and strain control using a fractional derivative inverse model[END_REF]. The model is constituted of two contributions. A quasi-static contribution related to the low-frequency behavior (f ≪ 1 Hz) and a dynamic contribution based on fractional operators and related to the frequency dependence of the hysteresis.

Quasi-Static Contribution

The easiest way to observe the quasi-static contribution is to plot the spontaneous polarization P , versus the electric field E, at very low frequencies (f ≪ 1Hz). At such frequency levels, we assume that Bloch wall movements will behave similarly to mechanical-like dry friction oscillations [START_REF] Ducharne | Low frequency modelling of hysteresis behaviour and dielectric permittivity in ferroelectric ceramics under electric field[END_REF][START_REF] Guyomar | Time fractional derivatives for voltage creep in ferroelectric materials: theory and experiment[END_REF]. Stressed by an external excitation, Bloch walls jump successively from one pinning defect to another until they reach a steady state where a minimum energy level is obtained. Each jump can be considered as a small mechanical dry friction, and consequently as a small amount of losses. A static (frequency-independent) equation based on its mechanical dry-friction counterpart has been established in order to consider this property. Basically, this equation describes how a major symmetrical hysteresis loop P(E) can be obtained by translating a hysteretic curve. The sign of the translation is equal to the sign of the time derivative of the polarization and the amplitude equal to the coercive field, Ec. In the time domain, the polarization can be written as:

P(t) = F E(t) -E c sgn dP(t) dt . (10.4)
Here, F(E) (or inversely F -1 (P )) represents the behavior of a nonlinear dielectric. Its mathematical description is as follows:

F(E(t)) = σ tan -1 E(t) γ . (10.5)
An experimental protocol has been defined to obtain this F(.) function. This procedure details are explained in detail in [START_REF] Ducharne | Low frequency modelling of hysteresis behaviour and dielectric permittivity in ferroelectric ceramics under electric field[END_REF]. The parameters γ , σ and the function F are obtained by fitting the tested samples' experimental anhysteretic curve to the analytical expression (Eq. 10.5). Unfortunately, due to its simplicity, Eq. 10.4 is limited to the description of symmetrical major hysteresis loops observed during steady state and for high-amplitude electrical field (E 0 ≫ E c ). Differences with experimental results will appear as soon as the excitation is no longer symmetrical (first polarization curve, minor hysteresis loops). Equation 10.4 is limited to just a single Bloch wall movement. To overcome this limitation and to get an accurate model independent of the experimental situation, we need to take into account a large number of domain wall dry frictions, as is the case experimentally. A distribution of mechanical dry frictions (called spectrum), characterized by their own coercive fields E ci and polarization P i in addition to their own weights, is taken into account to converge into much more precise simulation result:

P i = F E(t) -E ci sgn dP i (t) dt , (10.6)
k i=1 Spectrum(i) × P i (t) = P(t). (10.7)
Spectrum(i) represents the distribution of these basic dry frictions (wall movements). A deconvolution between the experimental first polarization curve and the F function allows to obtain this distribution. Large number of results as well as further information concerning the static model are available in [START_REF] Ducharne | Low frequency modelling of hysteresis behaviour and dielectric permittivity in ferroelectric ceramics under electric field[END_REF][START_REF] Guyomar | Time fractional derivatives for voltage creep in ferroelectric materials: theory and experiment[END_REF][START_REF] Zhang | Energy harvesting based on piezoelectric Ericsson cycles in a piezo ceramic material[END_REF][START_REF] Guyomar | High frequency bandwidth polarization and strain control using a fractional derivative inverse model[END_REF][START_REF] Rouleau | Application of Kramers-Kronig relations to timetemperature superposition for viscoelastic materials[END_REF].

Dynamic Contribution

Usually, frequency dependence in ferroelectric hysteresis models (Eqs. 10.6, 10.7) is considered by adding a viscous losses term to the quasi-static contribution. Such a consideration modifies Eq. 10.4 in the following way:

P(t) = F E(t) -E c sgn dP(t) dt + ρ dP(t) dt . (10.8)
Unfortunately, simple experimental test quickly demonstrates that this equation can only be used with a moderate accuracy on a restrained frequency bandwidth. If the comparison of simulation/measure for the amplitude-versus-frequency curve correctly fits the low-frequency range, the difference for the increasing frequency is important. The viscous losses term ρdP/dt, where ρ is a material constant, in Eq. 10.8 leads to an overestimation of the high-frequency component of the polarization signal. To overcome this problem, we need an operator that balances the low-frequency and the high-frequency component in a different way than a straight time derivative. Such operators exist in the framework of fractional calculus; there are so-called non-entire derivatives or fractional derivatives. The fractional derivative generalizes the concept of derivative to non-integer orders. The fractional derivative of a smooth (arbitrary) function g(t) is defined through a convolution between g(t) function and t α (t)/ Ŵ(1 -α) where Ŵ(α) is the gamma function; α, the order of fractional derivation; and (t) is the step Heaviside function (see [START_REF] Havriliak | A complex plane representation of dielectric and mechanical relaxation processes in some polymers[END_REF], while a summary of various representations can be found in [START_REF] Nguyen | Electrical and physical characterization of bulk ceramics and thick layers of barium titanate manufactured using nanopowders[END_REF]). After addition of fractional terms in our simulation, the model equations (Eqs. 10.6, 10.7) become:

P i = F E(t) -E ci sgn dP i (t) dt + ρ d α P i (t) dt α , (10.9)
k i=1
Spectrum(i) × P i (t) = P(t).

(10.10)

Model for the Complex Permittivity Under Weak Electric Field

At room temperature and as the electric field is weak, we admit that the real part ε(ω) and the imaginary part ε(ω) of the complex permittivity measured around an operating point are just frequency dependent. Both components of the complex dielectric permittivity are linked to each other by the well-known Kramers-Kronig relations [START_REF] Rouleau | Application of Kramers-Kronig relations to timetemperature superposition for viscoelastic materials[END_REF]. Debyes model has been one of the first models proposed to deal with the dielectric relaxation. But first experimental results have rapidly shown that corrections of that simple model will be necessary for a correct description. The Cole-Coles model introduces an improvement considering a distribution function for the relaxation time. This correction has been possible as a result of the extension of the Debyes model to a non-entire order. Cole-Coles model is indeed a frequency fractional dependence of the evolution of the complex permittivity. From the Cole-Coles model, ε(ω) and ε(ω) can be written as:

ε ′ (ω) = ε ∞ + Δε ′ 2 1 - sinh(βz) cosh(βz) + cos β π 2 , (10.11)
ε ′′ (ω) = Δε ′ 2 1 - sinh(β π 2 ) cosh(βz) + cos β π 2 , (10.12)
In the above expressions z = ln(ωτ ), Δε = ε s -ε ∞ , and β = (1 -α), where α shows the deformation of the semicircle arc in the Cole-Cole plot, i.e., it is the angle from the ε axis to the center of the semicircle arc. ε(ω) is related to the stored energy within the sample. ε(ω) is related to the dissipation (loss) of energy within the medium. According to the Cole-Coles model, the complex permittivity ε * may be written as:

ε * = ε ′ (ω) -iε ′′ (ω) = ε ∞ + Δε 1 + (iωτ ) α (10.13)
with Δε = ε s -ε ∞ and 0 <α= 1 -2σ/π ≤ 1, while σ is the angle with respect to the semicircle center. Various methods can be employed to determine this angle [START_REF] Nguyen | Electrical and physical characterization of bulk ceramics and thick layers of barium titanate manufactured using nanopowders[END_REF].

Relation Between High Excitation Amplitude Dielectric Hysteresis Fractional Model and the Complex Permittivity Cole-Cole Model

We start with a weak excitation contribution, sufficiently weak to assume that even for hysteretic material cases, there are linear relations between the dielectric polarization P and the electric field E. If the frequency is small and the polarization is harmonic:

P(t) = P 0 cos(ωt) (10.14)
the static (instantaneous) contribution of E is given by:

E stat (t) = P 0 cos(ωt) ε stat (10.15)
As the frequency is increasing, the dynamic contribution of the electric field is considered because of the fractional relation as it is the case in the hysteresis fractional model:

E dyn (t) = ρ d α P dt α . (10.16)
Here P is a harmonic-type waveform, so the dynamic contribution of the electric field can be expressed by its analytical solution:

E dyn (t) = ρP 0 ω α cos ωt + α π 2 .
(10.17)

Finally, for higher frequencies, the electrical field including all contributions is described by:

E(t) = ρP 0 ω α cos ωt + α π 2 + P 0 cos(ωt) ε stat (10.18)
In the complex representation, an electric field ℜ{E(t)}=E(t) and a corresponding dielectric coefficient ε can be written as follows:

E(t) = P 0 exp(iωt) 1 ε stat + ρω α exp iα π 2 , (10.19) ε = 1 1 ε stat + ρω α e iα π 2
, (10.20)

ε = ε stat 1 + ε stat ρω α i α .
(10.21)

Finally, α , where τ α = ρε stat .

ε = ε stat 1 + (iωτ )
(10.22)

If we assume ε(0) ≫ ε ∞ which is obviously our case as we work with piezoceramic where ε(0) ≈ 1800 and ε ∞ ≈ 100, the relation between Cole-Coles model and our high electric field amplitude model is clear:

ε = ε ∞ + ε(0) -ε ∞ 1 + (iωτ ) α) ≈ ε stat 1 + (iωτ ) α) .
(10.23)

Characterization and Experimental Validation

A soft PZT composition (P188 obtained from Quartz & Silice, France, Navy type II) has been tested in this study (Table 10.1). Cylindrical specimens (of 6.35 mm diameter and 4 mm height) are exposed to the external electric field. Constant temperature conditions and free mechanical properties of the sample are assumed.

High-Amplitude Electric Field Characterization

The excitation electric field is dimensioned to provide a maximal value close to 2 kV/mm. Under such level of excitation, we highly exceed the coercive field of the The polarization field is determined by integration of the charge measurement. Post processing, the numerical integration is performed. Figure 10.1 presents the experimental measuring setup especially developed for the high-amplitude electric field characterization. The role of the rod in the middle of the frame is to transmit the tested samples displacement to the displacement sensor. Even if the rod can have side effect because of its inertia, we have checked using laser vibrometer that the mechanical displacements are correctly transmitted by the system. For an electric field amplitude close to 2 kV, the frequency limitation of the high-voltage amplifier is close to 100 Hz for ceramics of 6 mm diameter. The excitation signal provided by the impedance analyzer consists of an AC sine of varying frequency with a maximal value corresponding to 1 V. The analyzer behaves as a perfect voltage source. The electrical current crossing the tested sample is monitored simultaneously. From the imposed voltage and from the measured current, the impedance analyzer provides instantaneously the paralleled voltageequivalent capacitance and the losses angle versus the frequency. The analyzer generates a linear sweep on a large frequency bandwidth available (40 Hz-40 MHz). In this study the maximum frequency tested has been reduced to 50 kHz to avoid the mechanical resonances which completely disturb the measure.

Experimental Results

Among others, an objective of this study is to check the accuracy of the fractional operators taken into account and the frequency dependence of the dielectric losses in a ferroelectric material. Both situations were tested then, the complex dielectric permittivity and the dielectric hysteresis dynamic under high level of external electric field.

As illustrated, in the model description in the above section, a theoretical frequency dependence relation can be established between Cole-Cole dielectric permittivity model and high excitation amplitude hysteresis model. In this experimental part, comparison of simulation/measure is exposed to validate first the viability of our theory and to check next if the dynamic parameters can be conserved through the scales. If this statement is confirmed, it means that whatever the scale, the dynamical behaviors are similar. It confirms too that a simple impedance meter analysis is enough to parameterize the materials ferroelectric dynamic behavior. Figure 10.2 shows comparisons of simulation/measure obtained under quasi-static conditions.

Figure 10.3 shows the quasi-static, high excitation amplitude, model parameters. Figure 10.4 shows similar comparisons for frequency varying from 2 MHz to 50 Hz. It is worth noticing that Fig. 10.4 exhibits good simulation/measure correlations on high-amplitude hysteresis loops. Unfortunately on such kind of measures, the highfrequency limit is rapidly reached. Beyond this maximal frequency, the current required to polarize the piezo ceramic is so high that the voltage amplifier is no longer suitable to provide it. One solution is to reduce the surface area of the tested sample, but even with a maximal reduction, we are rapidly limited. For better clarity, Fig. 10.5 illustrates the limitation of the first-order derivative model. In this figure, we have intentionally plotted only the high-frequency simulation results obtained for α = 0.53. However, the results obtained for f = 0.5 Hz and higher are similar to those obtained as α = 1. After validation of the high-amplitude excitation dynamic hysteresis model, it is time now to focus on complex dielectric model. Impedance meter gives the frequency evolution of the real and the imaginary part of the permittivity.

Cole-Cole plot can be obtained then and compared to those obtained with the Cole-Cole model. Cole-Cole model is configured using the same parameters to those giving good results with the hysteresis model. Unfortunately, the high-frequency limit of the impedance meter is relatively low. Indeed, as soon as the first mechanical resonant appears, the impedance of the tested sample changes radically, and the voltage source is no longer suitable to supply to the ceramic. The half circle Cole-Cole plot obtained is actually not even one third of a circle. The frequency bandwidth tested is sufficiently large to conclude that the comparisons obtained give very interesting results. The simulated curves obtained with the Cole-Cole model and the experimental ones fit almost perfectly and validate our expectations. The physical behaviors responsible for the dynamic dielectric losses are similar under high and weak electrical excitation field. They can then be modelled using the same operators. This observation is particularly interesting because it allows limiting the piezo ceramic characterization to the impedance analyzer (where all the dynamic model parameters can be set) and to anticipate the high electrical amplitude stress behavior in simulation. Impedance analyzer measures are quick and easy and can be repeated without special attention which is obviously not the case in high-amplitude hysteresis measure.

Conclusions

Fractional operators taken into account the dynamic dielectric losses through a piezo ceramic material under high electrical excitation were tested successfully in the past [START_REF] Guyomar | Time fractional derivatives for voltage creep in ferroelectric materials: theory and experiment[END_REF][START_REF] Weron | Probabilistic basis for the Cole-Cole relaxation law[END_REF][START_REF] Rekanos | Approximation of GrünwaldLetnikov fractional derivative for FDTD modeling of Cole-Cole media[END_REF][START_REF] Lewandowski | Fractional-order models: the case study of the supercapacitor capacitance measurement[END_REF]. This technique gives good accuracy on frequency dependence of hysteresis loops (see full black lines with α ≈ 0.5i nF i g .10.4). It also allows to envisage other dynamic manifestations such as creep behavior, ageing, and extended electrode contacts. In this article, the link between the fractional consideration used under high electric field amplitude in the limit of low frequency (required for hysteresis loop plots) and very low excitation level but high-frequency bandwidth measures was obtained using an impedance analyzer and discussed in the context of the Cole-Cole model. In both limits, it was confirmed that the dynamic losses exhibit the same physical origin, and consequently they could be modelled using the same operators and the same parameters. This notation is particularly interesting because it allows to reduce the piezo ceramic characterization to the impedance analyzer characterization. In this article, authors have focused their study to soft piezo ceramics. Future work will deal with other types of ferroelectric material such as polymers.

Fig. 10 . 1

 101 Fig. 10.1 Experimental measuring setup

Fig. 10 . 2 Fig. 10 . 3

 102103 Fig. 10.2 Comparison of simulation and measurement results in a quasi-static limit (f<50 mHz) for different initial polarizations (initial conditions for (a)-(c))

Fig. 10 . 4

 104 Fig. 10.4 Comparison of simulation and measurement results under dynamic conditions. Note that the convergence of the simulated (black line) for α ≈ 0.5 and experimental results (black points) and the lack of convergence (simulation results in red lines lower panels) for α = 1

Fig. 10 . 5

 105 Fig. 10.5 Comparison of simulation and measurement results in the Cole-Cole plot (ε(ω) versus (ε(ω))

Table 10 .

 10

	1 PZT (P188) material parameters
	Parameter	Symbol	Units	Typical value
	Density	ρ ′	10 3 kg m -3	7.7
	Poisson's constant	σ ′	-	0.3
	Curie temperature	T c	• C	340
	Dielectric permittivity	ε T 33 /ε 0	-	1850
	Piezoelectric coefficient	d 33	pC/N	425

Table 10 . 2 Table 10 . 2

 102102

	parameters	Simulation	Quasistatic parameters Dynamic parameters Symbol Va l u e Symbol Va l u e
			γ	1000	α 0.53
			σ	30	ρ	50,000

Acknowledgements

The authors would like to thank for the support from the Polish-French collaboration project (Polonium).


Luca Alasio
email: luca.alasio@sorbonne-universite.fr

Towards a new mathematical model of the visual cycle

come

Towards a new mathematical model of the visual cycle

Introduction and elements of the retinoid visual cycle

Photons reaching photoreceptors in the retina initiate the chemical process by which light is converted into electrical signals (phototransduction). A series of chemical transformations known as the visual cycle allows photoreceptors to return to the dark state and eventually re-iterate the process. This fundamental biochemical mechanism in rods and cones is activated by light reaching rhodopsin and causing retinals (vitamin A variants) to detach, as discussed in detail later. The visual cycle was first understood by George Wald, who later obtained the Nobel Prize in Medicine in 1967 for his discoveries.

The human retina is constituted by several types of cells arranged in multiple layers, however this project focuses on photoreceptors and the adjacent retinal pigment epithelium (RPE). Photoreceptor cells are the primary light-sensing cells in the retina. Photoreceptor proteins in the cells absorb photons, producing a change in the cell's trans-membrane potential. This triggers impulses that reach the brain through nerve layers connected to the optic nerve. Such process is referred to as phototransduction.

Rods, cones and their distribution

Two types of photoreceptors contribute directly to sight, namely rods and cones. Both rods and cones are specialised neuroepithelial cells capable of converting light into biological signals. Rods are narrower than the cones and distributed differently across the retina, which contains approximately 120 million rod cells and 6 million cone cells. Rods and cones share a similar structure: the most elongated parts, called rod/cone outer segments (ROS/COS), are in close contact with the RPE cell layer. [START_REF] Philip D Kiser | Chemistry of the retinoid (visual) cycle[END_REF] More precisely, rods consist of an outer segment (with respect to the center of the eye), a thin connecting cilium, an inner segment, cell body and synaptic region. The outer segment is composed of a stack of membrane discs surrounded by a separate plasma membrane [START_REF] Robert S Molday | Atp-binding cassette transporter abca4: molecular properties and role in vision and macular degeneration[END_REF]. The rod outer segment is roughly 25 µm in length and contained discs are roughly 16 nm thick and equally spaced [START_REF] William | Photoreceptor discs: built like ectosomes[END_REF] (see figure 1). The structure of cones is similar to that of rods, however their outer segment is shorter and its shape resembles a smoothed out cone. The chemical process for phototransduction is similar in rods and cones, however rods are extremely sensitive and can be triggered by a single photon, whereas cones are only sensitive to thousands of photons (e.g. in good light conditions). Cones require a larger number of photons to produce a signal. In humans, there are three different types of cone cells responding to different light wavelengths (blue, green and red). The centre of the retina is called macula and it is specialised for the central, high-resolution, colour vision in the presence of good light. Fovea and foveola are inner regions of the macula containing a high density of cones, spiking at the centre (see figure 2).

RPE cells

Retinal pigment epithelium (RPE) is a layer of cells supporting the neurosensory retina. An RPE cell adjacent to a photoreceptor provides it with nutrients, contributes to the visual cycle and phagocytosis of outer segments as part of the disc renewal process ( [START_REF] Robert S Molday | Atp-binding cassette transporter abca4: molecular properties and role in vision and macular degeneration[END_REF]), which takes approximately 10 days. Between 20 and 40 photoreceptors are adjacent to a single RPE cell. Such number varies with the retina since the distribution of rods and cones is not uniform and since the dimension of photoreceptors can vary slightly [START_REF] Philip D Kiser | Chemistry of the retinoid (visual) cycle[END_REF]. For example, the number of rods per RPE ranges from 29 in the periphery to 22 in the fovea [START_REF] Gao | Aging of the human retina. differential loss of neurons and retinal pigment epithelial cells[END_REF]. From our perspective, RPE cells have two crucial roles. Firstly, RPE guarantees photoreceptor cell renewal, which is essential because outer segments contain high levels of unsaturated lipids ( [START_REF] Nicolas G Bazan | Cell survival matters: docosahexaenoic acid signaling, neuroprotection and photoreceptors[END_REF]) subject to oxidation in the presence of light, (photo)reactive retinal, and high oxygen tension ( [START_REF] Clemens | Oxygen sensing in retinal health and disease[END_REF][START_REF] Caprara | From oxygen to erythropoietin: relevance of hypoxia for retinal development, health and disease[END_REF]). Perturbations of metabolite renewal and damaged component disposal can lead to exacerbated retinal degeneration ( [START_REF] Philip D Kiser | Chemistry of the retinoid (visual) cycle[END_REF]). Photoreceptors regenerate approximately 10% of their OS volume per day [START_REF] Brian | Phagocytosis of retinal rod and cone photoreceptors[END_REF]. The furthermost discs from the nucleus are shed, phagocytosed by adjacent RPE cells and replaced with newly formed outer segments. As a consequence, RPE cells receive a massive quantity of oxidised cellular debris ( [START_REF] Philip D Kiser | Chemistry of the retinoid (visual) cycle[END_REF]) which may not always be properly digested and hence tends to accumulate over several years. In particular, condensation compounds derived from all-trans-retinal (a photo-sensitive molecule similat to vitamin A) in RPE cells are toxic and cause oxidative stress ( [START_REF] Sparrow | The bisretinoids of retinal pigment epithelium[END_REF]). It has been proposed that cell migration or proliferation of peripheral RPE cells towards the macula can compensate for continuous loss of foveal RPE cells due to aging [START_REF] Priore | Age-related changes in human rpe cell density and apoptosis proportion in situ[END_REF].

Figure 1: On the left, an image of rod and cone cells (see [START_REF] Nelson | The Organization of the Retina and Visual System[END_REF]). On the right, two sections of the rod outer segment showing the structure of stacking membrane discs (see [START_REF] Gilliam | Three-dimensional architecture of the rod sensory cilium and its disruption in retinal neurodegeneration[END_REF]).

Figure 2: Distribution of rods and cones (left) and sub-regions of the center of the retina (right). Images adapted from [START_REF]Cmglee Zyxwv99[END_REF].

Visual cycle

The visual cycle is the process allowing a rod cells to return to the dark state after exposure to light (see figures 3 and 4). In particular, after photoexcitation, all-trans retinal dissociates from rhodopsin and is subsequently reduced to all-trans retinol on the cytoplasmic surface of disc membranes. All-trans retinol then reaches the adjacent RPE cells where it is converted to all-trans retinyl ester and subsequently to 11-cis retinol. 11-cis retinol is oxidized to 11-cis retinal it is transported back to the rod photoreceptor where it recombines with opsin in disc membranes to regenerate rhodopsin ( [START_REF] Robert S Molday | Atp-binding cassette transporter abca4: molecular properties and role in vision and macular degeneration[END_REF]).

The cone visual cycle differs from the rod visual cycle as it is faster and additional processes are involved ( [START_REF] Philip D Kiser | Chemistry of the retinoid (visual) cycle[END_REF]). After exposure to bright light, cones return to the previous state in 3-4 minutes, whereas rods take over 30 minutes ( [START_REF] Wang | The cone-specific visual cycle[END_REF]).

Specific molecules, proteins and enzymes

A large number of light sensitive molecules are needed to increase the probability of photon capture. These light-sensitive molecules must also couple to a receptor protein capable of initiating the required biochemical chain of events ( [START_REF] Philip D Kiser | Chemistry of the retinoid (visual) cycle[END_REF]). Let us list the main microscopic contributors to the visual cycle.

1. Vitamin A, in the form of all-trans-retinal, all-trans-retinol, all-trans-retinyl ester, 11-cis-retinol and 11-cis-retinal, plays a crucial role as a photosensitive molecule ( [START_REF] Philip D Kiser | Chemistry of the retinoid (visual) cycle[END_REF]).

2. Rhodopsin are membrane proteins, they are densely packed into specific regions of photoreceptors. Rhodopsin consists of two components, a protein molecule (an opsin) embedded in the lipid bilayer of cell membranes, and a photoreactive chromophore (retinal). The protein forms a pocket where retinal lies ( [START_REF] Daniel J Müller | Vertebrate membrane proteins: structure, function, and insights from biophysical approaches[END_REF]). There are approximately 10 8 rhodopsins per photoreceptor.

3. The enzymes responsible for the main steps of the visual cycle are RPE65 isomerohydrolases, 11-cisretinol dehydrogenases (11-cis-RDH), all-trans-retinol dehydrogenases (all-trans-RDH), lecithin retinol acyltransferases (LRATs). Their role will not be discussed further at this point (see [START_REF] Philip D Kiser | Chemistry of the retinoid (visual) cycle[END_REF] for further details). These three enzymes form a protein complex in RPE cells (see [START_REF] Strauss | The retinal pigment epithelium in visual function[END_REF])

4. Binding proteins and membrane transporters produced by photoreceptors allow hydrophobic retinals to diffuse in the aqueous phase and to cross membranes. Such proteins include interphotoreceptorbinding protein (IRBP), cellular retinaldehyde-binding protein (CRALBP), phosphatidylethanolamine (PE) and ATP-binding cassette (ABCA4). See [START_REF] Philip D Kiser | Chemistry of the retinoid (visual) cycle[END_REF] for further details.

As explaied in [START_REF] Wald | The molecular basis of visual excitation[END_REF][START_REF] Philip D Kiser | Chemistry of the retinoid (visual) cycle[END_REF], Wald discovered that the red color of rhodopsin was due to 11-cis-retinal, a compound derived from vitamin A. Moreover, photons cause a photochemical reaction transforming 11-cisretinal into all-trans-retinal. This is the first step in the sensing of light by the eye. Following photoactivation, rhodopsin decomposes into its protein and retinal components.

Different time and space scales.

The project considers phenomena taking place at significantly different space and time scales. This means that the derivation of a good model will be mathematically challenging. With regards to time scales, there are substantial differences:

• photobleaching of rhodopsins takes place in fractions of a second,

• the visual cycle takes several minutes (3-4 minutes for cones, 30 minutes for rods),

• regeneration of outer segments of rods takes approximately 10 days,

• macular atrophy progresses over decades.

The difference in spatial scales is equally dramatic:

• a membrane disc is about 10 nm thick,

• the rod outer segment is about 25 µm long,

• the macula has a diameter of approximately 5 mm.

Let us also note the large number of proteins and cells involved:

• number of rods in the retina: ranges between [1.1, 1.25] × 10 8 ,

• number of cones in the retina: 4.6 × 10 6 cone cells,

• number of discs per ROS (rod outer segment): approximately 1100,

• number of rhodopsins per rod photoreceptor cell: approximately 10 8 . Such figures are reported by BioNumbers [START_REF] Milo | Bionumbers-the database of key numbers in molecular and cell biology[END_REF].

Mathematical modelling: state of the art

In order to illustrate our mathematical approach, I propose a simplified model for the time and space evolution of the concentrations of all-trans-retinal and 11-cis-retinal in the visual cycle.

A first model

Consider a domain Ω ⊂ R 3 which is a perforated cylinder of radius R and height H, representing a rod outer segment. The "holes" in Ω represent the membrane discs, which is labeled by an index i = 1 . . . N discs ∼ 1000. They are supposed to be equi-spaced, centered around the axis of the cylinder, to have the same radius R disc < R and height H discs ∼ H 2 N -1 discs , so that their cumulative height is approximately H/2. Σ i indicates the boundary of each disc, Σ := Remark 2.1. The functions β j (j ∈ {0, 1, 2}) shall be determined in terms of reaction kinetics by the properties of the enzymes LRAT, RPE65 and 11-cis-RDH. As observed in [START_REF] Philip | Retinal pigment epithelium 65 kda protein (rpe65): An update[END_REF] and references therein, the efficiency of 11-cis-retinol production is dependent on the presence of retinoid-binding proteins suggesting that product release from the enzyme may in fact be rate-limiting. This suggests that the functions β i (•) may take the following Michaelis-Menten form:

β j (s) = V j max s K j + s . (1)
Our model, to be completed with suitable initial conditions, is the following:

• equations for the concentrations of all-trans-retinal u and 11-cis-retinal v in the rod outer segment Ω,

∂ t u -λ trans ∆u = 0 in Ω (2) ∇u • n = -α u K + u ≤ 0 on Γ (3) ∇u • n = ℓz rho ≥ 0 on Σ (4) ∂ t v -λ cis ∆v = 0 in Ω (5) ∇v • n = α v out K + v out ≥ 0 on Γ (6) ∇v • n = -vz ops ≤ 0 on Σ (7) (8)
• equations for the concentrations of rhodopsiins z rho and opsins z ops on each membrane disc Σ i ,

d dt z rho = -ℓz rho + vz ops on Σ (9)
d dt z ops = ℓz rho -vz ops on Σ (10) (11)
• equations for the concentrations of other retinoids in RPE cells (on the outer boundary Γ), Let us briefly analyse each equation in the model:

d dt u out = α u K + u -β 0 (u out ) on Γ (12) d dt u 1 out = β 0 (u out ) -β 1 (u 1 out ) on Γ ( 13
)
d dt u 2 out = β 1 (u 1 out ) -β 2 (u 2 out ) on Γ ( 14
)
d dt v out = β 2 (u 2 out ) -α v out K + v out on Γ (15)
• Equation ( 2) represents standard diffusion of all-trans-retinal in Ω, i.e. the idealised outer segment of a rod cell. The value of the diffusion coefficient depends on the properties of the cytoplasm and on those of all-trans-retinal and its binders.

• Equation (3) accounts for outflow of all-trans retinal through Γ. When all-strans retinal reaches the OS outer membrane it is transported out by IRBP. The limited carrying capacity is modeled by the logistic term.

• Equation ( 4) describes the inflow of all-trans retinal through Σ. In the presence of light, dissociation of retinal and opsin occurs, determining a source of all-trans-retinal in our equation. The constant ℓ vanishes in the absence of light.

• Equation ( 5) represents standard diffusion of 11-cis-retinal in Ω. It is analogous to equation ( 5).

• Equation ( 6) accounts for inflow of 11-cis-retinal through Γ. After various chemical tranformations in an adjacent RPE cell, 11-cis-retinal flows into the OS (Ω). IRBP acts once again as a transporter.

• Equation ( 7) describes the outflow of 11-cis-retinal through Σ. When reaching Σ, 11-cis-retinal combines with an opsin and form rhodopsin.

• Equation ( 9) represents the formation of rhodopsins in the presence of both opsins and 11-cis-retinal and their transformation into opsins in the presence of light.

• Equation (10) represents a reaction kinetics of first order for the variation of opsins. It is strictly related to equation ( 9) since sum of the concentrations of opsins and rhodopsins is constant over time.

• Equation ( 12) represents a reaction kinetics of first order for the variation of all-trans-retinol. It is associated to the flow in equation ( 6).

• Equation ( 13) represents a reaction kinetics of first order for the variation of all-trans-retinyl.

• Equation ( 14) represents a reaction kinetics of first order for the variation of 11-cis-retinol.

• Equation ( 15) is a first order reaction kinetics for the variation of 11-cis-retinal.It is associated to the flow in equation ( 6).

Remark 2.2. The sub-system consisting of equations ( 9) and ( 10), namely

d dt z rho = -ℓz rho + vz ops on Σ. d dt z ops = ℓz rho -vz ops on Σ.
can be simplified as follows:

z ops = C z -z rho , d dt z rho = -ℓz rho + v(C z -z rho ) = -(ℓ + v)z rho + C z v.
Furthermore, as a first order ODE for each point of Σ, it can be explicitly solved in terms of v:

z rho (t) = exp - t 0 (ℓ + v) dτ z rho (0) + C z t 0 v exp s 0 (ℓ + v) dτ ds . ( 16
)

Model improvement

There are several ways in which our model of the visual cycle can be improved. Here is an inexhaustive list:

1. All-trans-retinal is slightly hydrophobic and it can diffuse freely into the aqueous cytoplasm only in presence of a suitable binder. This aspect has not been considered in detail so far. On the other hand, all-trans-retinal can diffuse on the lipid membrane discs, which can be taken into account modifying our model.

2. The specific characteristics of the IRBP transporter have not yet been taken into account. IRBP is a large protein with two sites for all-trans-retinal and other two sites for 11-cis-retinal (preferably). This fact should lead to more complicated terms describing the flow inside and outside of Ω. A first step towards a better description involves the differentiation of rates of removal (namely the constant α) in equations ( 3) and ( 6), since the removal of 11-cis-retinal appear to be faster than that of all-transretinal, see [START_REF] Chen | Interphotoreceptor retinoid-binding protein removes all-trans-retinol and retinal from rod outer segments, preventing lipofuscin precursor formation[END_REF].

3. The chemical reactions involving retinals in the RPE cell shall be described in greater detail taking into account the action of enzymes such as RPE65, 11-cis-RDH, all-trans-RDH, LRATs. Additionally, the exchange of vitamin A isomers with the blood circulation must be taken into account in an improved model. Self-limiting effects in the chemical reactions are also possible (see e.g. [START_REF] Moiseyev | Rpe65 is the isomerohydrolase in the retinoid visual cycle[END_REF]).

4. In the limit of infinitely many membrane discs, an homogenisation problem arises: is it possible to derive an effective system of equations that captures the essential features of the visual cycle without keeping track or the single membrane disc? The powerful techniques presented in [START_REF] Gaudiello | Homogenization of the brush problem with a source term in l[END_REF] should be applied to get a mathematical homogenisation result in this direction.

5. The results in [START_REF] Td | Dark adaptation and the retinoid cycle of vision[END_REF] concerning dark adaptation should be reproduced and extended using our model.

6. In reality rhodopsin is not uniformly distributed on membrane discs, but it is rather arranged into densely packed granules in the discs, and absent on the rims. The model can be adapted to better capture such heterogeneity. To this end, we may redefine ℓ = ℓ ω (t, x) as a function of space, time and an event ω in a suitable probability space. This would allow to capture both spatial heterogeneity of discs and the stochastic nature of light irradiation.

Stationary solutions with constant illumination

As a first step in the analysis of system ( 2)-( 15), we consider the problem of boundedness and uniqueness of stationary solutions. In this subsection we assume that ℓ > 0 is given and constant in time and space. This corresponds to constant and homogeneous illumination of the photoreceptors for a very long time. Excluding time dependence, the system takes the following form:

λ trans ∆u = 0 in Ω, ( 17
) ∇u • n = -α u K + u on Γ, ( 18
) ∇u • n = ℓz rho on Σ, (19)
λ cis ∆v = 0 in Ω, ( 20
)
∇v • n = α v out K + v out on Γ, ( 21
)
∇v • n = -vz ops on Σ, (22)
z ops = C z -z rho on Σ, (23)
(ℓ + v)z rho = C z v on Σ, (24) 0
= α u K + u -β 0 u out on Γ, (25)
0 = β 0 u out -β 1 u 1 out on Γ, ( 26
) 0 = β 1 u 1 out -β 2 u 2 out on Γ, ( 27
) 0 = β 2 u 2 out -α v out K + v out on Γ. ( 28
)
Equations ( 23) and ( 24) provide explicit expressions for z rho and z ops in terms of v and given constants, in particular

z rho = C z v ℓ + v . ( 29
)
Equations ( 25)-( 28) imply v out = u. Keeping these two relations in mind, from equations ( 17)-( 22) we deduce that, for some constant

C v , u = C v -v. (30)
Altogether we have obtained the reduced system:

∆u = 0 in Ω, ( 31
)
∇u • n = -α u K + u on Γ, ( 32
)
∇u • n = ℓC z (C v -u) ℓ + C v -u on Σ. (33)
Notice that, in this simplified setting, system (2)-( 15) reduces to a single PDE defining a nonlinear boundary value problem for the Laplace equation.

Lemma 2.3. Problem (31)-( 33) admits at most one weak solution in the Sobolev space H 1 (Ω). Furthermore the solution satisfies the bilateral bound 0 ≤ u ≤ C v almost everywhere in Ω.

Proof. For convenience, we consider the following modified system:

∆u = 0 in Ω, ∇u • n = -α u K + |u| on Γ, ∇u • n = ℓC z (C v -u) ℓ + |C v -u| on Σ,
having the following weak formulation for any η ∈ H 1 (Ω):

Ω ∇u • ∇η dx = - Γ α u K + |u| η dσ + Σ ℓC z (C v -u) ℓ + |C v -u| η dσ. ( 34
)
Letting η = (u) -we easily see that any weak solution u is non-negative. Indeed we have the following retaion:

-

Ω |∇(u) -| 2 dx = Γ α (u) 2 - K + |u| dσ + Σ ℓC z (C v (u) -+ (u) 2 -) ℓ + |C v -u| dσ ≥ 0, (35)
and consequently (u) -= 0 almost everywhere in Ω. Similarly, letting η = (u -q) + (Stampacchia's truncation), we obtain

Ω |∇(u -q) + | 2 dx = - Γ α u(u -q) + K + |u| dσ + Σ ℓC z (C v -u)(u -q) + ℓ + |C v -u| η dσ. ( 36
) Choosing q = C v , it follows that Ω |∇(u -q) + | 2 dx + Γ α u(u -q) + K + |u| dσ + Σ ℓC z (u -q) 2 + ℓ + |C v -u| η dσ = 0, (37)
hence u ≤ C v almost everywhere in Ω.

Let us observe that

u K + |u| - ũ K + |ũ| (u -ũ) = K(u -ũ) + u|ũ| -|u|ũ (K + |u|)(K + |ũ|) (u -ũ) (38)
In particular

(u|ũ| -|u|ũ)(u -ũ) =      0 if uũ ≥ 0 -2uũ(u -ũ) if u > 0, ũ < 0 2uũ(u -ũ) if u < 0, ũ > 0 , (39)
which is always non-negative.

Analogously we have

C -u ℓ + |C -u| - C - ũ ℓ + |C -ũ| (u -ũ) = -ℓ(u -ũ) + (C -u)|C -ũ| -|C -u|(C -ũ) (ℓ + |C -u|)(ℓ + |C -ũ||) (u -ũ) (40)
In particular

((C -u)|C -ũ| -|C -u|(C -ũ))(u -ũ) =      0 if (C -u)(C -ũ) ≥ 0 -2(C -u)(C -ũ)(u -ũ) if C -u > 0, C -ũ < 0 2(C -u)(C -ũ)(u -ũ) if C -u < 0, C -ũ > 0 , (41)
which is always non-positive. We obtain uniqueness considering the difference of the weak formulations for two solutions u and ũ:

Ω ∇(u-ũ)•∇η dx = - Γ α u K + |u| - ũ K + |ũ| η dσ + Σ ℓC z C v -u ℓ + |C v -u| - C v - ũ ℓ + |C v -ũ| η dσ. ( 42
)
Choosing η = u -ũ, thanks to the previous computations, we have

Ω |∇(u -ũ)| 2 dx + Γ α K(u -ũ) 2 (K + |u|)(K + |ũ|) dσ + Σ ℓ 2 C z (u -ũ) 2 (ℓ + |C v -u|)(ℓ + |C v -ũ|) dσ ≤ 0, (43)
thus u = ũ almost everywhere in Ω. As a consequence of non-negativity, we also deduce that |u| = u, therefore the modified system coincides with the original one.

3 Modelling degenerative diseases, the case of Stargardt disease.

Stargardt disease (STGD) is the most common inherited macular degeneration, affecting approximately 1:10000 individuals, with an estimated annual incidence between 1.1 and 1.3 per million population (all ages) (see [START_REF] Cornish | The epidemiology of stargardt disease in the united kingdom[END_REF]). The disease is characterised by progressive macular atrophy due to accumulation of lipofuscin-like material; it affects both photoreceptors and retinal pigment epithelium. This leads to weakened central and detailed vision and decreased visual acuity. The disease (STGD) is recessively inherited, and it is caused by a mutation of the ABCA4 gene (the protein encoded by such gene is described in the next paragraphs).

The condition was first recognized in 1909 by Stargardt ([43]) and the term fundus flavimaculatus was later used by Franceschetti et al. ( [START_REF] François | Fundus flavimaculatus[END_REF]), who described "fishlike" yellow-white retinal lesions now called "flecks". Typical Stargardt disease with macular atrophy is usually diagnosed before the age of 20 with decreased visual acuity, weakened central and detailed vision that gradually worsens; fundus flavimaculatus is usually diagnosed later ( [START_REF] Sahel | Clinical characteristics and current therapies for inherited retinal degenerations[END_REF]). Diagnosis is based on fundus examination, full-field and multifocal or pattern ERG (electroretinogram), and fundus imaging such as fundus autofluorescence and OCT (Optical coherence tomography) ( [START_REF] Sahel | Clinical characteristics and current therapies for inherited retinal degenerations[END_REF]). Later ages of onset have been associated with a more favourable visual prognosis.

The visual cycle and the ABCA4 transporter

In rod cells, phototransduction is initiated when light activates rhodopsin (a light-sensitive receptor protein) in membrane discs by converting 11-cis retinal to its all-trans isomer. The ABCA4 (photoreceptor-specific ATP-binding cassette transporter 4) gene encodes for a protein that actively transports all-trans-retinal, a variant of vitamin A, across the disk membranes after photoexcitation of rhodopsin. The protein is localised at the rims of each disc of the outer segments of both rod and cone photoreceptors. Most of the all-trans retinal diffuses to the cytoplasmic surface of disc membranes where it is reduced to all-trans retinol and is converted to 11-cis retinal (via the visual cycle). However, a significant fraction of the all-trans retinal diffuses to the lumen side of the disc membrane where it reacts with PE (phosphatidylethanolamine) to form N-retinylidene-PE. This compound is trapped on the lumen side of the disc membrane. As discussed by Molday ( [START_REF] Robert S Molday | Atp-binding cassette transporter abca4: molecular properties and role in vision and macular degeneration[END_REF]), ABCA4 transports N-retinylidene-PE trapped on the lumen side of the disc membrane to the cyctoplasmic side consuming ATP.

Pathophysiology of STGD

According to Molday ( [START_REF] Robert S Molday | Atp-binding cassette transporter abca4: molecular properties and role in vision and macular degeneration[END_REF]), loss or diminished function of ABCA4 as an N-retinylidiene-PE transporter can explain most of the characteristic features exhibited by patients with Stargardt macular degeneration and abca4 knockout mice. Ineffective transport of N-retinylidene-PE from the lumen to the cytoplasmic side of the disc membrane determines an accumulation of N-retinylidene-PE and all-trans retinal in photoreceptor disc membranes, which can react to form the compound A2PE. It is also possible for all-trans retinal to re-associate with opsin and form a retinal opsin complex capable of triggering the visual cascade in a less efficient way ( [START_REF] Jäger | Opsin/all-trans-retinal complex activates transducin by different mechanisms than photolyzed rhodopsin[END_REF]). Photoreceptor outer segments are continually being renewed, indeed packets of aged disc membranes are phagocytised by adjacent RPE cells; the outer segment is completely renewed over a period of 10 days. Normally, the aged disc membranes are completely digested by the lysosomal system of the RPE cells. In this process, A2PE from outer segment discs can be transformed into A2E and phosphatidic acid, however RPE do not have enzymes capable of further degrading A2E. As a result, A2E and related side products progressively accumulate in RPE cells. The constitute the "flecks", i.e. fluorescent lipofuscin deposits observed in individuals with Stargardt diseases and in abca4 knockout mice. As explained by Molday [START_REF] Robert S Molday | Atp-binding cassette transporter abca4: molecular properties and role in vision and macular degeneration[END_REF], A2E is toxic for RPE cells since it can compromise subcellular organelles and an inhibitor of normal RPE degradative functions. Importantly, in the presence of oxygen and light A2E can be converted into free radical epoxides which can cause death of RPE cells. The death of RPE cells will result in photoreceptor degeneration and a loss in vision since RPE cells are crucial for photoreceptor cell survival. This results in lesions (atrophy) whose growth rates strongly depend on the size of the initial lesion ( [START_REF] Strauss | Progression of stargardt disease as determined by fundus autofluorescence in the retrospective progression of stargardt disease study (progstar report no. 9)[END_REF]). In a healthy human RPE, lipofuscin accumulates progressively throughout the first sixty years of life to occupy up to 19% of cytoplasmic volume; in later age it tends to plateau ( [START_REF] Michael E Boulton | Studying melanin and lipofuscin in rpe cell culture models[END_REF]). In Stargardt disease the RPE contains up to seven times more lipofuscin than normal ( [START_REF] Frank G Holz | Inhibition of lysosomal degradative functions in rpe cells by a retinoid component of lipofuscin[END_REF]).

Connections with other diseases

Mutations in ABCA4 are not only responsible for Stargardt disease but are also a major cause of cone-rod dystrophy (CORD), retinitis pigmentosa (RP), Best's disease ( [START_REF] Priore | Age-related changes in human rpe cell density and apoptosis proportion in situ[END_REF][START_REF] Frank G Holz | Inhibition of lysosomal degradative functions in rpe cells by a retinoid component of lipofuscin[END_REF]). The effects of the various mutations on ABCA4 protein function may be responsible for this phenotypic variability; nonsense or frame-shift mutations with protein truncation being more commonly involved in CORD or RP than in typical STGD ( [START_REF] Sun | Stargardt's abcr is localized to the disc membrane of retinal rod outer segments[END_REF]). Photo-oxidation products of A2E have been reported to activate the complement system suggesting that A2E and related compounds may serve as a trigger for age-related macular degeneration ( [START_REF] Zhou | Complement activation by photooxidation products of a2e, a lipofuscin constituent of the retinal pigment epithelium[END_REF]).

Modelling the visual cycle and STGD

A main goal of this project is to build a model for the evolution of the macular atrophy in Stargardt disease by means of adequate mathematical tools. Subsequently, the model may be perfectioned to support clinical results ans studies. The elementary interactions between different types of cells and proteins may be described by means of discrete or individual-based models, however a description at the continuous level (concentrations or population densities) is convenient since the number of single proteins and cells involved is of the order of hundreds of millions. In order to describe the evolution of the macular atrophy, cell populations corresponding to photoreceptors and RPE have to be modelled, hence one shall consider cell densities as functions of time and space. The photoreceptors and RPE cell densities will be obtained as a solution of an appropriate system of differential equations that will also include equations for the toxic compounds that accumulate in an unhealthy retina. In the first instance, the main components of our models can be described as follows:

• density of photoreceptors, described by a differential equation involving time, position in the retina, density of RPE cells near the photoreceptor position, concentration of all-trans retinal.

• density of RPE cells, described by a differential equation involving time, position, density of photoreceptors, concentration of all-trans retinal and lipofuscin.

• toxins, i.e. all-trans retinal and lipofuscin, both in photoreceptors and RPE. The concentrations of toxins will be described by differential equations that may include terms representing diffusion, transport, clustering (for the accumulation of lipofuscin), reaction and return to the vision cycle. The exact number of components needed to describe the action of toxins has to be discussed further.

• external factors such as light, nutrients and drugs might influence the cell activities and appear explicitly in the previous components of the system.

A very simplified mathematical model of the visual cycle has been proposed in [START_REF] Td | Dark adaptation and the retinoid cycle of vision[END_REF]: it includes evolution over time of the concentration of retinals, neglecting space dependence of the quantities involved.

4 Strategy and impact

Next modelling steps

I propose a road-map for the derivation and validation of a robust model from the visual cycle to STGD.

1. Understanding scales. A solid, qualitative understanding of all the relevant mechanisms involved in the visual cycle is necessary in order to start modelling a healthy retina. Physiological mechanisms of the (healthy) human retina have been studied in the literature and their knowledge is well-established.

It is desirable to have the same level of information for an unhealthy retina, however this is not always available since the bio-molecular study of macular degeneration is rather recent. Like other complex phenomena arising in Life Sciences, the retina displays a great variety of mechanisms at very different scales. It is crucial to identify which phenomena matter the most at a given scale. The most significant mechanisms will be included in the models. I plan to focus on the following aspects, with the support of specialists of the retina:

• protein scale, typical dimensions: nanometres, minutes; I will collect quantitative information about the behaviour of proteins of the vision cycle and enzymes, interaction with membranes, active transport, diffusion within and across cells (see [START_REF] Chen | Formation of all-trans retinol after visual pigment bleaching in mouse photoreceptors[END_REF]),

• cell scale, typical dimensions: micrometres, days; I will concentrate on lipofuscin accumulation, phagocytosis, clustering mechanism [START_REF] Höhn | Lipofuscin: formation, effects and role of macroautophagy[END_REF][START_REF] Nociari | Lipofuscin accumulation into and clearance from retinal pigment epithelium lysosomes: Physiopathology and emerging therapeutics. Lysosomes: Associated Diseases and Methods to Study Their Function[END_REF][START_REF] Family | Dynamic scaling of lipofuscin deposition in aging cells[END_REF], toxicity and triggering of apoptosis [START_REF] Suter | Age-related macular degeneration: the lipofuscin componentn-retinyl-n-retinylidene ethanolamine detaches proapoptotic proteins from mitochondria and induces apoptosis in mammalian retinal pigment epithelial cells[END_REF][START_REF] Wenzel | Molecular mechanisms of light-induced photoreceptor apoptosis and neuroprotection for retinal degeneration[END_REF],

• tissue scale, typical dimensions: millimetres, years; I will study the interaction of RPE and photoreceptors as layers of cells (see [START_REF] Strauss | The retinal pigment epithelium in visual function[END_REF]), focusing on space inhomogeneity possibly leading to free boundaries in the model.

2. Model derivation. Once processes at different scales are well understood, I shall proceed and merge all the information in order to build a first, non-specialised model. This preliminary model might present several hybrid features, such as a combination of discrete and continuous components accounting for different mechanisms from the molecular level that of retinal layers. The preliminary model will be further refined as it shall become as simple as possible without losing any substantial information. A subsequent step consists in specialising the preliminary model into two models describing the healthy and the unhealthy retina. The specialised models are likely to take the form of a system of ordinary and partial differential equations, whose unknowns will be the density of photoreceptor cells, the density of RPE cells and the concentration of important chemicals such as A2E and lipofuscin in each type of cells. Nutrients and light will be suitably accounted for by means of additional terms in the equations. Finally, the model will be well-suited to include the effects of possible treatments, which can be included adapting the precise form of the equations and adding ad hoc parameters if necessary. For example, proliferation and migration of RPE or stem cells from a therapy might be introduced in the form of source and diffusion terms respectively, and the effect of drugs competing to replace the toxins might result in an additional equation in the system. In summary, the key steps are the following:

• derivation of the non-specialised model and related system of differential equations,

• refinement to a specialised model for the healthy retina,

• refinement to a specialised model for the unhealthy retina,

• addition of external factors (e.g. light, nutrients, stem cells and drugs).

3. Analysis and approximation of the models. Once the derivation is complete, I will proceed to analyse the models using the tools of Mathematical and Numerical Analysis. The role of mathematics is to understand the analytical structure of the new models, to explain the qualitative behaviour of solutions and to describe their fundamental properties rigorously and precisely. Nonlinear effects are to be expected, hence both core techniques of the analysis of differential equations (such as fixed-point methods, a priori estimates) and problem-specific tools will be very useful in the analysis. Two elements will make the analysis of the models particularly challenging are, firstly, the presence of clustering effects (concerning lipofuscin) and nonlinearities, possibly leading to integro-differential equations, secondly the need to study the long-time asymptotic behaviour of the system (relevant for predictions), which is often a difficult problem. I will also work on the numerical analysis of the model, finding methods to approximate and solve the system efficiently. The goal is to obtain an approximation scheme that is stable, consistent, and convergent and this can be achieved, for example, by means of finite difference and finite element methods. Preliminary numerical simulations will be performed at this stage. Uncertainty in the parameters and data of the model are expected, thus I will do a sensitivity analysis in order to identify the most significant factors affecting the uncertainty of the model predictions, as well as uncertainty quantification, dealing both with epistemic uncertainties (due to incomplete knowledge of the retina) and uncertainties (due to imperfect measurements). Finally, I will validate the model against a dataset for the progression of STGD (not used in the derivation) by estimating and tuning the parameters in order to improve the predictive performance: different regimes might be determined by dissimilar rates of progression and severity of the disease. Other analytical and numerical challenges are to be expected and will depend on the exact form of the equations. The main steps will include the following:

• well-posedness study (existence, uniqueness, regularity and continuous dependence),

• study of long-time asymptotics and stationary states,

• study of numerical approximation schemes and their convergence,

• preliminary simulations,

• sensitivity analysis and uncertainty quantification,

• validation, parameter estimation and tuning.

4. Interaction with ophthalmologists and medical applications. An in silico analysis through simulations will provide new insights in the evolution of STGD and other related diseases and the results will be discussed with ophthalmologists in order to get feedback and make the most effective use of such new tools. I hope to achieve significant extension of the current knowledge and to targeted simulations designed to answer specific questions of medical relevance. Finally, I will focus on the simulation of proposed therapies, such as drugs reducing the negative effects of toxins, regeneration of RPE thanks to stem cells or gene therapy. In short, the last part of the project will focus on:

• optimised numerical simulations,

• targeted simulations and knowledge extension,

• study and simulation of proposed therapies, customised simulations.

Research methodology and approach

Models can be designed in order to describe a phenomenon at different levels of detail, however simplifications are necessary given the great number of cells and molecules involved. In a first step the model should capture the typical qualitative behaviour of STGD, and subsequently move on to study the quantitative behaviour. The derivation and study of the model(s) will proceed by stepwise refinement, gradually improving the effectiveness of the description simultaneously reducing redundancy of details. As a starting point, it is possible to consider an individual-based model for toxins, their production and accumulation in a photoreceptor-RPE cell pair. This preliminary model should combine deterministic and stochastic features (see also [START_REF] Family | Dynamic scaling of lipofuscin deposition in aging cells[END_REF] for the Statistical Physics viewpoint): the work of David J. Aldous [START_REF] David | Deterministic and stochastic models for coalescence (aggregation and coagulation): a review of the mean-field theory for probabilists[END_REF] on models for coalescence will be useful in this context. Given the extremely large number of cells and molecules involved, individual based models are not ideal, whereas a continuous model (in space and dime) taking into account generation and transformation rates of the density of all-trans retinal and the related compounds is more likely to produce better results. Besides coagulation, diffusion and transport terms might also play a decisive role in capturing all the features of the visual cycle and STGD. In the passage from individual-based models to continuous models I will adapt the techniques developed by Juan J. L. Velázquez and collaborators [START_REF] Bodnar | An integro-differential equation arising as a limit of individual cell-based models[END_REF]. Finally, compartmental models originating from epidemiology and population dynamics might be used to describe cells passing from a state to another (e.g. healthy/unhealthy/apoptotic), as well as proteins clustering into various larger compounds which may be described as separate populations.

From the point of view of Mathematical Analysis, modelling accumulation of lipofuscin by means of coagulation-fragmentation equations represents a challenging novelty both from the theoretical and the applied perspective. The interaction of different components will affect the overall behaviour of the system. Previous work in this direction includes the results of Herbert Amann [START_REF] Amann | Coagulation-fragmentation processes[END_REF] and Marie Doumic [START_REF] Doumic | Eigenelements of a general aggregation-fragmentation model[END_REF], which are conveniently oriented towards long-time asymptotics. With respect to Numerical Analysis, methods developed by Yvon Maday (reduced-basis) and Albert Cohen (wavelet methods) might play an important role in our numerical schemes.

Earlier works concerning the development of mathematical models for the retina and its physiology are illustrated in [23]. It should be noted that the work in [START_REF] Yr Malashenko | Mathematical model of substance active transport via cytoplasmic membrane of microbial cell along concentration gradient[END_REF][START_REF] Paul A Roberts | Mathematical and computational models of the retina in health, development and disease[END_REF][START_REF] Paul A Roberts | Predictive mathematical models for the spread and treatment of hyperoxiainduced photoreceptor degeneration in retinitis pigmentosa[END_REF][START_REF] Vega | Anomalous angiogenesis in retina[END_REF][START_REF] Hoyle | Generative mathematical modelling to demonstrate virtual simulations of neovascular age related macular degeneration[END_REF][START_REF] Sala | Mathematical modelling and simulation of ocular blood flows and their interactions[END_REF] has led to models for different retinal mechanisms such as active transport via cytoplasmic membrane, retinitis pigmentosa, anomalous angiogenesis, neovascular age related macular degeneration and ocular blood flows.

Impact, originality and interdisciplinary collaboration

This project will promote the interaction of the communities of mathematicians and ophthalmologists, with the intent of developing innovative tools for the study of STGD. The idea of a mathematical model capable of predicting STGD progression is completely new. Only recently the disease has been studied on large groups of patients in vivo and over a long period of time (e.g. in the study ProgStar). Furthermore, quantitative bio-molecular and genetic studies are still ongoing and currently new variants of the disease are being discovered. A robust mathematical model capable of predicting the evolution over time of the area affected by macular atrophy is highly desirable and it will contribute to the study of multiple degenerative retinal diseases which require a deeper understanding. In a long term perspective, mathematical models and simulations provide powerful tools to study the disease progression and to put new therapies to the test in silico. The project is intrinsically interdisciplinary since a combination mathematical, biological and medical insights are necessary at all stages. Besides the mathematical model itself, the project outcomes will include codes for simulation, new directions for disease study, management and therapy. Throughout this project I expect myself and the other researchers involved to gain expertise thanks to joint study and results, achieving a deeper understanding of the other disciplines and techniques involved, which will establish a new connection between Mathematics and Ophtalmology.

Figure 3 :

 3 Figure3: A schematic representation of a rod outer segment illustrating the presence of rhodopsins/opsins on membrane discs and their role during phototransduction in combination with all-trans-retinal and 11-cisretinal (image adapted from[START_REF] Krishnavedala | The chemical reactions involved in a visual cycle[END_REF]).

Figure 4 :

 4 Figure 4: A schematic representation of the different steps in the visual cycle (image adapted from [4]).

Figure 5 :

 5 Figure 5: Organization of rod photoreceptors in the mammalian retina (see [13]). (A) Vertical section of a vertebrate retina stained with toluidine blue. The retina contains five classes of neurons (photoreceptors, bipolar cells, horizontal cells, amacrine cells, ganglion cells), and is divided into three nuclear layers (outer nuclear layer (ONL), inner nuclear layer (INL), ganglion cell layer (GCL)) and two synaptic layers (outer plexiform layer (OPL), inner plexiform layer (IPL)); (B) Schematic of a vertebrate rod photoreceptor consisting of the long, light-sensitive outer segment (OS) enclosed by the retinal pigment epithelium (RPE) which is connected via a small intracellular bridge, the connecting cilium (CC), to the metabolically active inner segment (IS) including the basal body complex (BBC; shaded in blue) and the apical region. (C-H) Transmission electron micrograph of a photoreceptor-connecting cilium; (C) The CC connects the OS and IS. The BBC consists of the basal body (BB) and its centriole (Ce); (D) Cross-section of the CC showing the characteristic (9 × 2 + 0) microtubule array; (E) RPE with pigment granules and the Bruch's membrane (BM); (F) OS composed of stacks of tightly packed membrane discs; (G) Higher magnification view of the ciliary apparatus; (H) Synaptic ribbon (SR) with invaginations of postsynaptic horizontal cells (HC) and bipolar cells (BC). Scale bars: 100 µm (A), 200 nm (D,F,H), 500 nm (E,G).

N

  discs i=1 Σ i , and Γ indicates the external boundary of the cylinder. The exterior normal vector is denoted by n (regardless of the different components of the boundary, by slight abuse of notation).The following functions and parameters are considered:u concentration of all-trans-retinal in Ω v concentration of11-cis-retinal in Ω z ops concentration of opsins on Σ z rho concentration of rhodopsins on Σ u out concentration of all-trans-retinol that has left Ω through Γ (towards RPE) u 1 out concentration of all-trans-retinyl ester in RPE u 2 out concentration of 11-cis-retinol ester in RPE v out concentration of 11-cis-retinal flowing into Ω ℓ rate of rhodopsin activation / all-trans-retinal production α rate of transport of retinal by IRBP β 0 function describing production of all-trans-retinyl ester in RPE β 1 function describing production of 11-cis-retinol in RPE β 2 function describing production of 11-cis-retinal in RPE λ trans diffusivity of all-trans-retinal in Ω λ cis diffusivity of 11-cis-retinal in Ω K Michaelis constant

Figure 6 :

 6 Figure 6: Schematic representation of a vertical section of Ω. The boundaries Γ and Σ are coloured in blue and black respectively.

Acknowledgements I would like to thank Prof. B. Perthame for his constant support throughout this project and I would like to thank Prof. I. Audo, Dr. M. Nassisi and Dr. E. Nandrot for several useful discussions. This project 15 has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie grant agreement No 945298-ParisRegionFP, supported by the Paris Region Fellowship Programme.