Camille Krewcun

Laurent Sarry

Nicolas Combaret

Émilie Péry

Fast simulation of stent deployment with plastic beam elements

Coronary stent deployment is a reference cardiology intervention, used to treat atherosclerosis and prevent heart attacks. The outcomes of the intervention highly depend on the accuracy of the stent apposition, which could benefit from per-operative prediction tools. In this paper, we propose a fast and mechanically realistic 3D simulation of a coronary stent expansion. Our simulation relies on the finite element method and involves serially linked beam elements to model the slender geometry of a stent. The elements are implemented with a non-linear elasto-plastic behavior, describing realistically the complex deformation of a balloon-expandable stent. As a proof of concept, we simulated the free expansion of a coronary stent. The simulation output was compared with micro-CT data, acquired experimentally during the device expansion. Results show that the plastic beam model is able to reproduce successfully the final geometry of the stent. In addition, the use of 1D elements allows to achieve a significantly lower computational time than for equivalent literature simulations, based on 3D elements. This preliminary work highlights the compatibility of our method with clinical routine in terms of execution time. Further developments include the application of the method to more advanced simulation scenarios, with the addition of a personalized artery model.

I. INTRODUCTION

The coronary artery disease results from a physiological aging process causing a progressive narrowing (stenosis) of the artery lumen. In the case of coronary arteries, symptoms include an 'acute' phase resulting in myocardial infarction (i.e. heart attack). In this study, we focus on the Percutaneous Transluminal Coronary Angioplasty (PTCA) procedure with stent deployment. The intervention consists in reducing the stenosis by inflating a balloon and expanding a metallic endoprosthesis (stent). When the balloon is deflated, the stent keeps its deformed shape and stays permanently inside the vessel, to support the artery wall. The outcomes of this intervention highly depend on the stent geometry and on the accuracy of the deployment. Consequently, numerical simulation is a valuable prediction tool to minimize the probability of complications.

We find in literature an important number of studies using the Finite Element Method (FEM) to simulate the deployment of coronary stents. A commonly mentioned objective is to study the influence of the stent geometry over the global mechanical behavior, 1 Camille Krewcun, Laurent Sarry, and Émilie Péry are with the Institut Pascal, Université Clermont Auvergne, CNRS, SIGMA Clermont, F-63000 Clermont-Ferrand, France. Email: camille.krewcun@uca.fr 2 Nicolas Combaret is with the institut Pascal, CHU Clermont-Ferrand Université Clermont Auvergne, CNRS, SIGMA Clermont, F-63000 Clermont-Ferrand, France. in order to optimize the stent design. Example of such works were published in [START_REF] Migliavacca | Mechanical behavior of coronary stents investigated through the finite element method[END_REF] and [START_REF] Petrini | Numerical investigation of the intravascular coronary stent flexibility[END_REF], in which the free (unconstrained) deployment of a stent is simulated with different geometries. A free expansion is also simulated in [START_REF] Beule | Realistic finite element-based stent design: The impact of balloon folding[END_REF] in order to study the influence of the balloon model on the global simulation. Studies such as [START_REF] Gijsen | Simulation of stent deployment in a realistic human coronary artery[END_REF] and [START_REF] Zahedmanesh | Determination of the influence of stent strut thickness using the finite element method: implications for vascular injury and in-stent restenosis[END_REF] additionally introduced an artery model in the simulation, allowing to assess the impact of the stent geometry on the artery wall stress. As priority is given to a high precision, almost all the simulated objects in the aforementioned studies (stent, balloon, artery) are discretized finely with volumetric finite elements. Although this allows to optimize the stent design prior to its production, the corresponding models can not be used directly in a clinical environment, for instance to perform patient-specific simulation. Indeed, the computational times required by 3D element discretization, and standard commercial codes are prohibitive compared to the duration of the corresponding clinical interventions (for instance, an execution time of 48 hours is mentioned in [START_REF] Beule | Realistic finite element-based stent design: The impact of balloon folding[END_REF]).

In this paper, we propose a faster simulation routine to model accurately the deployment of a balloon-expandable stent. Our objective is to reach an execution time compatible with clinical routine, so that our simulation method could be used in practice as a relevant prediction tool (typically in patient-specific applications). We decide to base our model on the discretization of the stent geometry by 1D serially linked beam elements. Similar work was proposed by Čanić, Tambača et al. in [START_REF] Tambača | Mathematical modeling of vascular stents[END_REF] and [START_REF] Čanić | Cardiovascular stents as pde nets: 1d vs. 3d[END_REF], where the stent structure is modeled by an assembly of 1D rod elements. A limitation of these studies is that the proposed mechanical behavior is only elastic, therefore corresponding to auto-expandable stents. Nevertheless, the computational time mention in [START_REF] Tambača | Mathematical modeling of vascular stents[END_REF] is significantly faster (with an order of magnitude of 1 s). We decide to use beam elements in a similar way, in order to model the more complex (nonlinear) deformation undergone by a balloon-expandable stent. The use of connected beam elements has already been proved efficient in medical simulation to model slender structures, such as catheters for navigation in [START_REF] Duriez | New approaches to catheter navigation for interventional radiology simulation[END_REF], endovascular coils for brain aneurysm in [START_REF] Dequidt | Interactive simulation of embolization coils: Modeling and experimental validation[END_REF], flexible needles in [START_REF] Duriez | Interactive simulation of flexible needle insertions based on constraint models[END_REF] and [START_REF] Adagolodjo | Inverse real-time finite element simulation for robotic control of flexible needle insertion in deformable tissues[END_REF], or vascular networks [START_REF] Peterlík | Modeling and real-time simulation of a vascularized liver tissue[END_REF]. The use of beam elements in these studies allows to achieve low computational times (up to interactive simulations), but once again is limited to elastic deformations.

II. METHOD

A. Serially-linked beam model

The first step of our method consists in discretizing the stent structure, by considering its wire-like geometry as a series of straight slender segments. In order to take advantage of this particular shape, we model the stent struts using serially connected beam elements. Beams are discretized into 2 nodes, each described by 6 Degrees Of Freedom (DOFs): three degrees for the node's position, and three for its orientation. Linking beam elements allows to refine the discretization of the curvilinear stent parts. Typically, the number of beam elements discretizing a given strut length increases with curvature. A limitation of this model is that it is primarily intended to represent simple structures, in which different beam sequences are not connected. In the case of vascular stents, we need to connect more than 3 beam elements at the same node, while conserving the information of each beam's orientation. This particular issue is addressed using Lagrangian constraints, as developed in section II-F.

B. Mesh conception

Although the various commercialized coronary stents represent a wide range of shapes, all can be decomposed into basic elements. In this way, coronary stents are characterized by two elementary features:

• stent rings are constituted by the periodic repetition of a single geometrical pattern. The dimensions of this particular pattern are sufficient to reconstruct the ring. • rings are unified in the global stent structure by connectors, which are repeated periodically along the stent length between each pair of rings.

Consequently, only the dimensions of the ring base pattern and the connectors are required to reconstruct the entire stent structure. In order to retrieve these geometric features, we perform a series of X-ray acquisitions on a crimped device.

Considering the dimensions of coronary stents (and vascular stents in general), we use micro-CT imaging to obtain a sufficiently high resolution.

With the stent key dimensions, we can reproduce the rings' base pattern and connectors. From these elements, we can take advantage of the aforementioned stent periodicity to simply replicate circumferentially the base pattern, and axially the resulting ring and the connectors. Concretely, the conception is carried out using a Computer Aided Design (CAD) software (Solidworks TM ). All structures are designed using edges, which can then be converted to straight beam elements. As the association of a beam element to an edge requires to compute the orientation of the edge's nodes, we only use CAD to design a flat version of the global stent mesh. This flat version is then wrapped into a cylindrical shape when the beam elements are created.

C. Finite Element Method implementation

Our implementation of the FEM routine relies entirely on the open-source simulation framework SOFA1 , which was first mentionned in a publication in [START_REF] Allard | Sofa -an open source framework for medical simulation[END_REF]. The fundamental principle of the simulation workflow in SOFA is to rely on Newton's second law of motion, which can be simply written as:

Ma t+h = f t+h (x, v), (1)
where M is the mass matrix of the discretized mechanical system (i.e. if the system is described by n DOFs, then M is a n × n matrix), x is the global position vector, v the velocity, and a the acceleration. Subscript t indicates that the variable is considered at time t, and h is the time step.

Adopting the Backward Euler implicit formulation introduced in [START_REF] Baraff | Large steps in cloth simulation[END_REF] we obtain the following system:

(M -hB -h 2 K) A ∆v t+h = hf t (x, v) + h 2 Kv t b , (2)
where K = ∂f ∂x is the element stiffness matrix and B = ∂f ∂v the damping matrix. A more detailed version of the global simulation process can be found in [START_REF] Faure | SOFA: A Multi-Model Framework for Interactive Physical Simulation[END_REF]. The last linear system (A∆v t+h = b) is solved in ∆v t+h at each time step. From ∆v t+h and the Backward Euler scheme, we then retrieve x t+h making the system progress forward in time.

In [START_REF] Petrini | Numerical investigation of the intravascular coronary stent flexibility[END_REF], we use non-specific expressions for M and B: namely a uniform mass model (i.e. M is diagonal and all its elements are equal), and a Rayleigh damping model (i.e. B = r M M + r K K, r M and r K representing the Rayleigh mass and Rayleigh damping respectively). On the contrary, the expression of the stiffness matrix K highly depends on the finite element discretization and mechanical behavior. An important part of this study consists in proposing a novel approach to compute K, using a nonlinear elasto-plastic beam model.

D. Elasto-plastic behavior

A critical issue in our study is to handle correctly the complex deformation undergone by balloon expandable stents. During balloon inflation, endoprostheses are deformed from a crimped rest configuration to comply to the deploying balloon membrane. When the balloon is deflated, internal stresses resulting from the deformation do not tend to bring the stent back to its rest shape. This phenomenon indicates a highly plastic deformation, which must be appropriately simulated.

In the following, we describe briefly the main principles in the implementation of a plastic beam model, by focusing on the computation of the associated stiffness matrix K. For more details on the application of the FEM to beam theory, or in the definition of K, we refer the reader to [START_REF] Przemieniecki | Theory of Matrix Structural Analysis[END_REF] or [START_REF] Bonnet | Analyse des solides défomables par la méthode des éléments finis[END_REF].

Linear elasticity

In the case of linear elasticity, the stiffness matrix of a single beam element is computed as follows. First, we define as Ω(x, y, z) the volume represented by the beam element. Coordinates x, y, and z represent material coordinates, so that

x ∈ [0, l], y ∈ [-h 2 , h 2 ] and z ∈ [-w 2 , w 2 ]
respectively, where l represents the beam length, and h and w the height and width of its section. We use the approximation that h and w remain constant during the deformation. We note u(x, y, z, t) the continuous displacement field in the element. As mentioned in II-A, beams are discretized into 12 DOFs. Thus, we can also define the nodal displacement vector of a beam as ũ = (w x,i , w y,i , w z,i , θ x,i , θ y,i , θ z,i ), where i ∈ 1, 2 indicates which of the beam nodes the DOFs are referring to. The relation between u and ũ is described by shape functions, defining the beam continuous medium by interpolation of the nodal DOFs. In this study, we use 3D Timoshenko beam functions, as expressed by Bazoune and Khulief in [START_REF] Bazoune | Shape functions of three dimensional timoshenko beam element[END_REF]. The Timoshenko beam model allows notably to take into account shear deformations and to describe more accurately beams for which the section dimensions are comparable to the length [START_REF] Bazoune | Shape functions of three dimensional timoshenko beam element[END_REF]. In our case, this will be the case in the stent parts presenting a high curvature. The shape functions can be summarized in the form of a 3 × 12 matrix N giving:

u(x, y, z, t) = N(x, y, z)ũ(t). (3)
The expression of N is explicitly given in [START_REF] Bazoune | Shape functions of three dimensional timoshenko beam element[END_REF].

From the displacement field, we can compute the strain tensor using the small strain formulation. We represent the second order strain tensor as a 6 × 1 vector, using Voigt notation.

(x, y, z, t) = Su(x, y, z, t),

with

S =           ∂ ∂x 0 0 0 ∂ ∂y 0 0 0 ∂ ∂z 0 1 2 ∂ ∂z 1 2 ∂ ∂y 1 2 ∂ ∂z 0 1 2 ∂ ∂x 1 2 ∂ ∂y 1 2 ∂ ∂x 0           and =         x,x y,y z,z y,z x,z x,y         . (5)
Combining ( 3) and ( 4), we obtain:

(x, y, z, t) = B e (x, y, z)ũ(t), (6)
where B e = SN is a 6 × 12 matrix. Using Hooke's law for isotropic linear elasticity, we can relate the stress tensor σ to the strain:

σ(x, y, z, t) = C (x, y, z, t). (7)
It can be noted that in [START_REF] Čanić | Cardiovascular stents as pde nets: 1d vs. 3d[END_REF], σ and are second order tensors, imposing C to be expressed as a fourth order tensor. Using again Voigt notation, and σ can be expressed as 6dimensional vectors, while C is reduced to a 6 × 6 matrix. The corresponding expression for C can be found in [START_REF] Przemieniecki | Theory of Matrix Structural Analysis[END_REF]. From the Principle of Virtual Work, we can express the beam stiffness matrix K in integral form, using the quantities above, as:

K = Ω B e T CB e dΩ. (8)
The resulting K is a 12 × 12 matrix. Detail of the computation from the continuous expression of the Principal of Virtual Work can be found in [START_REF]Theory of Matrix Structural Analysis[END_REF]. We stress out the fact that matrices B e and C are constant over time, meaning that K could be computed once and for all at the beginning of the simulation. As elements of B e are at most second degree polynomials, we could also compute analytically the integral. This is namely used in [START_REF] Duriez | New approaches to catheter navigation for interventional radiology simulation[END_REF], [START_REF] Dequidt | Interactive simulation of embolization coils: Modeling and experimental validation[END_REF], and [START_REF] Peterlík | Modeling and real-time simulation of a vascularized liver tissue[END_REF] to save computational time. Nevertheless this becomes irrelevant when replacing the linear elasticity law by a nonlinear plastic behavior. Consequently, we concretely compute K by using Gaussian reduced integration over each beam element.

Plasticity

In plasticity theory, the constitutive law (i.e. the stressstrain relation) can be decomposed into two phases:

• an elastic phase, during which the internal stress depends linearly on the strain, as in ( 7), • a plastic phase, during which part of the strain can be distinguished as resulting from plastic energy dissipation, and the stress-strain relation becomes nonlinear. The transition from the elastic to the plastic phase is described by a yield criterion f : R 6 → R, defined on the stress space. In this study, we use the Von Mises yield function:

f (σ) = σ e (σ) -σ 0 , (9)
where σ e (σ) = (σ e1 + σ e2 )

1 2 , (10)
with

σ e1 = 1 2 (σ x -σ y ) 2 + (σ y -σ z ) 2 + (σ z -σ x ) 2 , (11)
and

σ e2 = 3 σ 2 xy + σ 2 yz + σ 2 zx . (12)
In ( 9), σ 0 is a material constant referred to as the yield stress. In a 1D tensile test, it corresponds to the stress value above which the material starts to deform plastically. It is determined experimentally for a given material.

The yield criterion is used as follows:

• if f (σ) < 0, the stress state σ corresponds to an elastic deformation, • if f (σ) > 0, σ is actually an unacceptable stress state, i.e. the deformation can not be considered as elastic anymore. A new realistic stress state must be computed. The computation of an acceptable plastic stress state relies on the decomposition of the strain into a plastic and an elastic component:

= el + pl . (13)
In the following, we adopt a model of perfect plasticity, as described in [START_REF] Hughes | Theoretical foundation for large-scale computations of nonlinear material behavior[END_REF]. The term 'perfect' refers to the fact that the yield stress σ 0 remains constant during the plastic phase. Consequently, the yield surface described by f (σ) = 0 in the stress space also remains unchanged.

The evolution of the plastic strain is described by a flow rule, of the form:

d pl = dλ ∂g ∂σ . (14)
The scalar value dλ is known as the plastic multiplier, and g : R 6 → R is a plastic potential defined on the stress space.

In our case, we use an associative flow rule, meaning that the plastic potential is considered equal to the yield criterion:

g(σ) = f (σ). ( 15
)
This approximation is generally acceptable for metals, but gives incoherent results with more complex elasto-plastic materials (such as soil, or concrete).

As the constitutive law becomes nonlinear in plastic phase, it is not possible to use [START_REF] Čanić | Cardiovascular stents as pde nets: 1d vs. 3d[END_REF] to compute the new stress state at each time step anymore. Additionally, the expression of K given in [START_REF] Duriez | New approaches to catheter navigation for interventional radiology simulation[END_REF], relying on the fact that C is constant, is no longer relevant. An alternative is to integrate the differential elasto-plastic stress-strain relation, in order to compute the new stress state incrementally, at each time step. A solution for this well known issue is to adopt an integration scheme based on the implicit return algorithm. Although the idea was introduced before, we refer the reader to [START_REF] Hughes | Theoretical foundation for large-scale computations of nonlinear material behavior[END_REF] for a detailed implementation, and to Simo and Taylor's study [START_REF] Simo | Consistent tangent operators for rateindependent elastoplasticity[END_REF] for a thorough review of this procedure. We give a brief description of the method below.

We want to compute a new stress increment ∆σ so that σ t+h = σ t + ∆σ. The computation is decomposed into two steps. At first, we compute an incremental elastic predictor, as:

∆σ el = C∆ . (16)
Depending on the sign of f (σ t + ∆σ el ), we consider two scenarios:

• if f (σ t + ∆σ el ) < 0, the new stress state σ t+h actually results from a purely elastic deformation, and the elastic predictor is accurate. Therefore ∆σ = ∆σ el . This scenario is also chosen if ∇ T f (σ t )∆σ el < 0 (unloading case). • if f (σ t + ∆σ el ) > 0, the elastic predictor does not provide an acceptable stress state. We must couple it with a corrective term, in order to obtain a plastic stress state. As described in [START_REF] Hughes | Theoretical foundation for large-scale computations of nonlinear material behavior[END_REF], the orthogonal return algorithm consists in searching for a corrective stress increment which, summed to the elastic predictor, gives the final stress state σ t+h . The term 'return' refers to the fact that the direction of this corrective stress increment is given by the unit normal to the yield surface. When using the Von Mises criterion, the regularity of the function implies that this direction is the radial direction. This is known as the radial return algorithm.

As given in [START_REF] Hughes | Theoretical foundation for large-scale computations of nonlinear material behavior[END_REF], we obtain the following update for the new stress state:

σ t+h = σ 0 σ e (σ el t+h ) s el t+h + Tr(σ el t+h )I, (17)
where σ el t+h = σ t + ∆σ el and s el t+h is the corresponding deviatoric stress (s el t+h = σ el t+h -Tr(σ el t+h )I). Following the consistency condition (∇f T dσ = 0), the plastic multiplier increment is computed as:

∆λ = ∇f T σ el t+h ∇f T C∇f . ( 18
)
Once the stress state σ t+h is known, we may finally compute the resulting internal forces, and linearized stiffness matrix, so that they are taken into account in the global mechanical system (2). The new internal forces are simply computed by integration:

f int (σ t+h ) = Ω B e T σ t+h dΩ. (19)
The linearized version of the stiffness matrix (or tangent stiffness matrix) K t can be expressed in a similar analysis way as for ( 8), leading to:

K t = Ω B e T dσ d B e dΩ. ( 20
)
Starting from Hooke's law, we have:

dσ = C d el = C(d -d pl ) dσ = C(d -dλ ∂g ∂σ ).
Through calculus, and using the consistency condition:

∇f T dσ = 0, (21)
we obtain an equivalent differential relation between stress and strain in plastic deformation:

dσ = C ep d , (22)
with

C ep =   C - C ∂f ∂σ ∂f ∂σ T C ∂f ∂σ T C ∂f ∂σ   . ( 23
)
This allows to explicitly express the tangent stiffness as:

K t = Ω B e T C ep B e dΩ. (24)

E. Expansion forces

In order to simulate the balloon expansion, we apply radial forces on each mesh node. These forces are incremented linearly in time, up to a maximal intensity value f i where i represents the index of a node. We approximate this intensity by considering that the inflation pressure of the balloon is directly transmitted to the stent inner surface S int . The global force intensity F applied on the stent can then be related to the inflation pressure P by:

F = P × S int = P × w × h × l tot , ( 25
)
where l tot stands for the mesh total length, computed by summing the lengths of all the elements. Then, we associate to each node a scalar value α i corresponding to the ratio of the lengths of the beam elements it belongs to, over 2l tot (as each beam length is counted twice: once for each of its nodes). α i gives an estimation of which proportion of the stent length is affected by node i.

We simply use these coefficients to distribute global force intensity over the nodes, leading to:

f i = α i F. ( 26
)
In order to reproduce the balloon inflation and deflation, we increment and decrement linearly the intensity of the nodal forces, over a given number of time steps. The direction of theses forces remains constant (along radial direction) during the deformation.

F. Constraints

As mentioned in the beginning of this paper, we need to ensure nodal cohesion between three or more beam elements, in certain parts of the mesh (for instance at stent ring connections). This is achieved by making use of Lagrangian constrained dynamics, as introduced in [START_REF] Witkin | Physically based modeling: Constrained dynamics[END_REF]. This method was specifically described on beam elements in SOFA by Sorokin et al. in [START_REF] Sorokin | Filogen: A modelbased generator of synthetic 3d time-lapse sequences of single motile cells with growing and branching filopodia[END_REF], to simulate filopodia. Once again, we use the same model without modifications, and refer the reader to the aforementioned article for more details. Briefly, the addition of constraints involves an additional step in the numerical solution;

• At first, (2) is solved without being constrained, leading to a new position (free motion). • In a second time, a modified constrained system is solved:

A∆v t+h = b + hH T λ, (27)
where H is a Jacobian containing the constraint directions, and λ contains the Lagrange multipliers associated to each constraint. In SOFA, this system is typically solved iteratively using a Gauss-Seidel algorithm.

III. RESULTS

As a proof of concept, the methodology presented above was entirely carried out on a Presillion Plus (Cordis, Fremont, California, U.S.) coronary stent.

A. Mesh conception

A first run of micro-CT images was acquired on the crimped stent, to retrieve its geometrical base features. The acquisition was performed with a micro-CT scanner (Explore CT 120 R , GE Healthcare, Chicago, Illinois, U.S.) with a predefined routine of 700 projections. The dimensions of the stent base patterns were measured on the resulting 3D image, and a corresponding flat mesh was conceived accordingly, as shown in Fig. 1. The studied stent can be decomposed into two sinusoidal ring patterns, used alternatively over a total of 13 rings. The global structure was discretized with 1956 beam elements (2005 nodes), and its cohesion assured with 85 bilateral constraints.

B. Simulation parameters

The flat mesh was imported into SOFA, and wrapped into a cylindrical shape. The radial expansion forces applied to each node were computed for an inflation pressure of P =3.0 MPa (about 30 atm), and the FEM routine was run with ∆t =1 µs. The simulation was launched on a common workstation (Intel Xeon E3-1270 v5 (8CPUs) 3.6 GHz, 16 GB RAM, Windows10 64 bits), without any form of optimization (single thread CPU execution). The execution time for the whole simulation was around 1 hour.

C. Experimental validation

In order to assess the precision of our proof-of-concept model, we performed three other micro-CT acquisitions, on a Presillion Plus deployed with an inflation pressure of 12, 16 and 20 atm, respectively. Stent expansion was unconstrained (i.e. no artery phantom was used), and the parameters were the same as for the crimped stent acquisition. Quantitative evaluation was achieved by comparing specific geometrical features on two pointsets, extracted from the simulation output and from the micro-CT 3D reconstruction, respectively. The simulation pointset simply corresponds to the mesh nodes, while the micro-CT pointset was extracted by converting the volumetric representation of the stent to a binary skeleton. We performed rigid registration between the two pointsets. In a first step, rigid transform was carried out to superimpose the principal axes of inertia and centers of gravity. In a second step, we used a rigid rotation around the common axis to minimize the Euclidean distance between simulated points and their closest neighbors in the CT pointset. 3D rendering of the registration output is given in Fig. 3 for a crimped and a deployed configurations.

The first available feature is the Euclidean distance between each pair of coupled points. We measured an average distance of 0.17 mm, which is approximately 10 times smaller than the macroscopic dimensions of the stent (the diameter and length's order of magnitude being typically 1 mm). In the same way, we computed the radial distance between each pair of points. We obtained an average value of 80 µm, which is the same order of magnitude as the strut thickness (the stent section dimensions being (50 × 50) µm). The last geometrical feature is the deployment radius, computed for each point as the radial distance with the stent rotation axis. Fig. 2 shows the mean radius in function of the axial position, along the stent centerline. We observe that the expansion radius is accurately estimated by the simulation around 0, with a minimum difference of 2.29 × 10 -3 mm between the two models, for an axial position of -0.28 mm. On the contrary, the radius is underestimated at the stent extremities, with a maximum difference of 0.210 mm for an axial position of -4.7 mm. We note that the free deployment of the stent leads to an heterogeneous radius along the centerline, which is correctly reproduced by the simulation. In practice, when a coronary stent is deployed inside an artery, its final radius is imposed by the artery wall.

IV. DISCUSSION

The first point we want to mention is that the model we propose for balloon-expandable stents is able to reproduce the free expansion of a stent, in a reasonable amount of time. The comparison of the simulation output with experimental micro-CT data, based on basic geometrical features, showed that the deployed configuration could be accurately reproduced in terms of radial expansion. In addition, the achieved execution time is significantly smaller than for similar simulations using volumetric elements and standard commercial codes. Among the studies mentioned in introduction and also modeling stent free expansions, but with 3D elements ([1], [START_REF] Petrini | Numerical investigation of the intravascular coronary stent flexibility[END_REF] and [START_REF] Beule | Realistic finite element-based stent design: The impact of balloon folding[END_REF]), only De Beule et al. in [START_REF] Beule | Realistic finite element-based stent design: The impact of balloon folding[END_REF] mention an execution time: less than 48 hours, for simulations run with ABAQUS large deformation analysis on a in-house-built high performance computing cluster. Although their discretization is particularly dense, we can expect similar computational times for [START_REF] Migliavacca | Mechanical behavior of coronary stents investigated through the finite element method[END_REF] and [START_REF] Petrini | Numerical investigation of the intravascular coronary stent flexibility[END_REF] in which the simulation also relies on ABAQUS large deformation analysis. Our 1D beam model is likely to be 10 times faster (48 times faster than [START_REF] Beule | Realistic finite element-based stent design: The impact of balloon folding[END_REF]), while running on a regular workstation (single CPU thread). One perspective of our work is to run an optimized version of our code on a more powerful computer (Intel Xeon E5-2630 v4 2.20 GHz (20 cores)). As several tasks can be run in parallel over the beam elements (namely the computation of the internal forces, the tangent stiffness matrix, and plastic stress states), we are confident in being able to divide the execution time by at least 10.

Nevertheless, the results also exhibit non-negligible limitations. We can see in the bottom picture of Fig. 3 that the stent length in the final configuration is shorter in the simulation than in the micro-CT imaging (9.35 mm against 11.1 mm). Even with similar expansion diameters, this unexpected shortening creates an axial shift in the whole structure. This phenomenon could be explained by various approximations of our model, the most likely of which are listed below. The perfectly plastic model is an obvious approximation for the plastic behavior. Even if hardening phenomena are not crucial considering the unique loading/unloading cycle undergone by the stent during its deployment, we could straightforwardly extend our model by taking into account the evolution of the yield criterion, as proposed in [START_REF] Hughes | Theoretical foundation for large-scale computations of nonlinear material behavior[END_REF]. Under the same expansion forces, such a modification should notably stabilize the stent deformation in comparison with perfect plasticity. The balloon simulation is another potential issue. As the balloon inflation and deflation are only simulated through incremental forces, we do not render directly the equilibrium that, in reality, takes place between the balloon elastic membrane and the coronary stent. Typically the implementation of this equilibrium could limit the shortening phenomenon mentioned above. The last clear limitation regarding our results is the fact that the default implementation of the Newton-Raphson algorithm in SOFA (used to solve (2)) restricts the procedure to 1 iteration per time step. Although this does not represent an issue when dealing with linear elastic behaviors, the introduction of a nonlinear plastic behavior should benefit from a fully iterative procedure. Although we restricted the accuracy assessment of this proofof-concept study to basic geometrical criteria, future work should ensure that the precision of our model is compatible with clinical use.

V. CONCLUSION

In this paper, we present a complete methodology to realistically simulate the plastic deformation of a balloonexpandable stent, in a reasonable time. The efficiency of the algorithm relies on an original discretization of the sent by 1D serially linked beam elements. Even if this study is focused on coronary stents, we highlight the fact that -to the best of our knowledge -this geometrical representation could be applied to other models of (uncovered) vascular stents, for instance designed to treat the carotid, femoral, or iliac arteries. Additionally, we experimentally assessed the output of the simulation using basic geometrical criteria, thus establishing the feasibility of the method and the compatibility of our model with clinical routine. As this work remains preliminary, it is meant to be extended towards more realistic applications. Consequently, the objectives of our future investigations will be the enhancement of the present beam model and the introduction of a constraining artery model in the simulation. Using a personalized model (for instance retrieved from Optical Coherence Tomography endovascular imaging), the simulation could be integrated in a clinical routine, in the form of a reliable visualization tool.

Fig. 1 .

 1 Fig. 1. Flat mesh reproducing the Presillion Plus coronary stent structure. The mesh was conceived with CAD, from micro-CT-acquired geometrical data.

Fig. 2 .

 2 Fig.2. Comparison of the Presillion Plus expansion radius between the simulation output (orange curve) and a micro-CT ground-truth acquired with P = 12 atm (blue curve). The result is given in function of the stent axial position, along its rotation axis.

Fig. 3 .

 3 Fig. 3. 3D rendering of the rigid registration between the simulation output and the micro-CT acquisition. The CT volumetric reconstruction is displayed plain and in grey, while the simulation mesh nodes are superimposed as red points. The top image corresponds to the crimped configuration (P = 0 atm), and the bottom image to a deployed configuration (P = 12 atm in the micro-CT acquisition).

www.sofa-framework.org

ACKNOWLEDGMENT

The authors would like to thank Arnaud Briat and Leslie Mazuel from the Imagerie Moléculaire et Stratégies Théranostiques (IMOST) laboratory (UMR 1240 Inserm / UCA), for their assistance in the acquisition of the micro-CT images.


Pierre-Olivier Lamare

Nina Aguillon

Jacques Sainte-Marie

Jérôme Grenier

Hubert Bonnefond

Olivier Bernard

Gradient-based optimization of a rotating algal biofilm process

Keywords: Adjoint-based optimization, PDE Control, microalgae, biofilm, biofuel, Rotating Algal Biofilm

Microalgae are microorganisms which have been only recently used for biotechnological applications, especially in the perspective of biofuel production. Here we focus on the shape optimization and optimal control of an innovative process where the microalgae are fixed on a support. They are thus successively exposed to light and dark conditions. The resulting growth can be represented by a dynamical system describing the denaturation of key proteins due to an excess of light. A Partial Differential Equations (PDE) model of the Rotating Algal Biofilm (RAB) is then proposed, representing local microalgal growth submitted to the time varying light. An adjoint-based gradient method is proposed to identify the optimal (constant) process folding and the (time varying) velocity of the biofilm. When applied to a realistic case, the optimization points out a particular configuration which significantly increases the productivity compared to a base case where the biofilm is fixed.

Introduction

Algae are microorganisms which were so far rarely used for biotechnological applications. In the last decade, their potential for production of bioproducts has been highlighted for addressing different markets such as pharmaceutics, food, feed, or even, at longer horizon, green chemistry and biofuels [START_REF] Wijffels | An outlook on microalgal biofuels[END_REF]. These organisms are generally cultivated under a planktonic form in open or closed photobioreactors. But the cell concentration remains low since for higher biomass light cannot pene-This paper was not presented at any IFAC meeting. Email addresses: pierre-olivier.lamare@mines-paristech.fr (Pierre-Olivier Lamare), aguillon@ljll.math.upmc.fr (Nina Aguillon), jacques.sainte-marie@inria.fr (Jacques Sainte-Marie), jerome.grenier@student.ecp.fr (Jérôme Grenier), hubert.bonnefond@inalve-conception.com (Hubert Bonnefond), olivier.bernard@inria.fr (Olivier Bernard).

trate through the very turbid medium [START_REF] Posten | Microalgae Biotechnology[END_REF]. Cell biomass thus only represent a tiny fraction of the liquid, generally around 0.1%. A high energy demand is then required by the harvesting and dewatering phases, and deeply limit the economic and environmental sustainability of these classical algal culturing process [START_REF] Lardon | Life-cycle assessment of biodiesel production from microalgae[END_REF][START_REF] Mata | Sustainability considerations of biodiesel based on supply chain analysis[END_REF]. In the recent years, surface attached microalgae biofilms have emerged [START_REF] Gross | Development of a rotating algal biofilm growth system for attached microalgae growth with in situ biomass harvest[END_REF]. They consist in growing microalgal cells fixed on a support and forming a biofilm. A biofilm is an assemblage of microbial cells that are irreversibly associated with a surface and enclosed in a matrix of extracellular polymeric substances. Biofilms are ubiquitous in nature, covering all kinds of surfaces, such as rocks, plants and sediments in seawater and freshwater environments. Algae in a biofilm is much more concentrated and straightforwardly harvested by a simple scrapping mechanism. Recently, innovative processes have been proposed where the algal biofilm is fixed on a rotating support so that cells are successively submitted to illuminated and dark periods (see Figure 1). The main advantage of this so called Rotating Algal Biofilm (RAB) is that it can better manage the way light is supplied to the microalgae and the subsequent dynamics of the photosystems harvesting light can be better exploited. When exposed during a too long period to high light intensity, some key proteins in the photosystems are denatured and must be resynthesized, which penalize the photosynthetic activity [START_REF] Molina Grima | A study on simultaneous photolimitation and photoinhibition in dense microalgal cultures taking into account incident and averaged irradiances[END_REF]. The purpose of this work is to propose a first model of this new process, accounting for the dynamical character of the photosynthetic apparatus when submitted to varying light. The microalgae growth rate and the related productivity are given by the Han model representing the dynamics of the photosystems (Section 2). Roughly speaking productivity accounts for how much the biomass is growing. Then, a Partial Differential Equations (PDE) model based on the former derivations is designed for the RAB (Section 3). This nonlinear infinite dimensional model is further used to derive a gradient based strategy to estimate, depending on the Photon Flux Density (PFD) the optimal belt folding ratio (ratio of belt which is enlightened) and velocity (Section 4) with respect to a productivity criterion. This latter method is powerful for large system with a large number of control variables and sufficiently general to be applied to a wide variety of control optimization problem with PDEs constraints.

Originally, the adjoint method has been developed within the optimal control community with the seminal work of Lions [START_REF] Lions | Optimal Control of Systems Governed by Partial Differential Equations[END_REF]. Therefore, various control problems are tackled with this method as ramp metering [START_REF] Reilly | Adjoint-based optimization on a network of discretized scalar conservation laws with applications to coordinated ramp metering[END_REF], air-traffic control [START_REF] Bayen | Adjointbased control of a new Eulerian network model of air traffic flow[END_REF], inverse pressure design and drag minimization problems [START_REF] Nadarajah | A comparison of the continuous and discrete adjoint approach to automatic aerodynamic optimization[END_REF]. Very recently, the method has been applied to the state and parameter estimation for traffic flow and overland flow system in [START_REF] Nguyen | State and parameter estimation in 1-D hyperbolic PDEs based on an adjoint method[END_REF]. Finally, there exist alternative approaches for gradient-based optimization of biological systems such as [START_REF] Passino | Biomimicry of bacterial foraging for distributed optimization and control[END_REF].

Here we use the classical methodology of gradient-based optimization but in the non-trivial framework of a very complex, nonlinear, realistic, and innovative model. The closest works in terms of the control problem are [START_REF] Bayen | Adjointbased control of a new Eulerian network model of air traffic flow[END_REF][START_REF] Work | Convex formulations of air traffic flow optimization problems[END_REF]. In [START_REF] Work | Convex formulations of air traffic flow optimization problems[END_REF] the velocities vary in space and time, so that a relaxation technique can be used to transform the problem into a convex problem in the control variable. However, in the RAB the velocity cannot be space dependent due to the the rigidity of the conveyor belt. Finally, the main difference of our work compared to [START_REF] Bayen | Adjointbased control of a new Eulerian network model of air traffic flow[END_REF], beyond the application, is the treatment of the control constraints. In [START_REF] Bayen | Adjointbased control of a new Eulerian network model of air traffic flow[END_REF], the cost function is augmented to penalize the violation of the constraints. Here, we favoured a projection step during the descent algorithm to satisfy the constraints.

Photoinhibition Dynamics

Photosystems are the functional units associated with the light part of photosynthesis. They are made of protein complexes and absorb light to generate a flux of energy and electrons. There are two types of photosystems in series PSI and PSII (see [START_REF] Han | Photosynthesis-irradiance response at physiological level: A mechanistic model[END_REF][START_REF] Wu | A model integrating fluid dynamics in photosynthesis and photoinhibition processes[END_REF]). PSII is generally assumed to play the main role in the photosynthetic dynamics. In the Han model, PSII is described by three different states: open and ready to harvest a photon (A), closed while processing the absorbed photon energy (B), or inhibited if two many photons have been absorbed simultaneously (C). Their respective evolutions satisfy the following dynamical system

dA(t) dt = -σI(t)A(t) + B(t) τ (1)
dB(t) dt = σI(t)A(t) - B(t) τ + k r C(t) -k d σI(t)B(t) (2) dC(t) dt = -k r C(t) + k d σIB(t) , (3)
where A, B, and C are the relative frequencies of the three possible states and I is a continuous time-varying light signal In the sequel we drop the t argument to alleviate the equations. Besides, the sum of A, B, and C is equal to 1 allowing to rewrite system (1)-(3) as

dA dt = -σI(t) + 1 τ A + (1 -C) τ (4)
dC dt = k d - k r k d + σI(t) C + σI(t) (1 -A) . (5)
We complete system (4) and ( 5) with initial conditions

(A(0, •), C(0, •)) = A 0 , C 0 ∈ E , (6)
where

E = (A, C) ∈ R 2 + |A + C ∈ [0, 1] .
Photosynthetic production is described by the transition between open state and closed state. Excitation is assumed to occur at a rate σI, with σ [m 2 µE -1 ] being the effective cross section of the photosynthetic unit (PSU), whereas deexcitation is assumed to occur at a rate of 1 τ , with τ [s] the turnover time of the electron transport chain. Photoinhibition occurring at high light irradiance corresponds to the transition from closed state to inhibited state. This process is assumed to occur at a rate of k d σI, with k d [-] a damage constant. The reverse transition from inhibited state to closed state accounts for the repair of damaged photosynthetic units by enzymatic processes in the cell, a mechanism that is assumed to occur at a constant rate k r [s -1 ]. The photosynthesis rate is thus proportional to σIA.

The dynamics of resting state A is fast compared to the dynamics of the photoinhibition C (due to the presence of the small multiplicative parameter k d in [START_REF] Gross | Development of a rotating algal biofilm growth system for attached microalgae growth with in situ biomass harvest[END_REF], see Table 1 below for an example of a possible value for this parameter). Since we focus on light variations at time scales larger than second, we can do a slow-fast approximation, and assume that A rapidly reaches its steady state [START_REF] Hartmann | The effect of photosynthesis time scales on microalgae productivity[END_REF] and focus on the dynamics of the photoinhibition rate C. To lead this analysis we use the singular perturbation theory (see e.g. [START_REF] Khalil | Nonlinear Systems[END_REF]). We define

K 1 = ε -1 k r (7) K 2 = ε -1 k d , (8)
where ε is a small parameter of magnitude order of k r and k d . Using parameters ( 7) and ( 8) equation ( 5) reads

dC dt = εK 2 - K 1 K 2 + σI(t) + σI(t) (1 -A) . (9)
Then, with the change of time

t = εt , (10)
and by defining A t = A(t), C t = C(t), and I t = I(t) we get the following system

ε d A d t = -σ Ĩ + 1 τ A + 1 -C τ (11) d C d t = K 2 - K 1 K 2 + σ Ĩ C + σ Ĩ 1 -A , (12)
with the initial conditions

A(0) = A(0) = A 0 (13)
C(0) = C(0) = C 0 . ( 14
)
Setting ε = 0 in [START_REF] Nguyen | State and parameter estimation in 1-D hyperbolic PDEs based on an adjoint method[END_REF] we get

A = 1 -C τ σ Ĩ + 1 . ( 15
)
Then, we substitute A in (12) by its expression given in [START_REF] Wu | A model integrating fluid dynamics in photosynthesis and photoinhibition processes[END_REF] we get the reduced dynamics for

C d C d t = K 2   -    K 1 K 2 + τ σ Ĩ 2 τ σ Ĩ + 1    C + τ σ Ĩ 2 τ σ Ĩ + 1    . ( 16
)
Let us define the boundary layer correction as

A (t) = A t - 1-C( t) τ σ Ĩ(t)+1 . We have d A(t) d t = d A t d t - d d t 1 -C t τ σ Ĩ t + 1 .
Therefore, we have

d A(t) dt = -σ Ĩ t + 1 τ A t + 1 -C t τ -ετ σ d Ĩ t d t τ σ Ĩ t + 1 -2 + ε τ σ Ĩ t + 1 d C t d t . Using A t = A(t) + 1-C( t) τ σ Ĩ( t)+1
and setting ε = 0 which freezes t at 0, we get the boundary layer model as

d A(t) dt = -σ Ĩ(0) + 1 τ A(t) . (17)
We have the following result.

Lemma 1 Let A and C be the solution to (4)-( 6), then for any fixed T > 0 there exists ε 0 such that 0 < ε < ε 0 and

C (t) -C (t) = O (ε) (18)
A (t) - 1 -C (t) τ σI(t) + 1 -A (t/ε) = O (ε) (19)
hold uniformly for t ∈ [0, T ], where C denotes the solution of equation

d C dt = -α(I) C + β(I) , (20)
where

α(I) = k d τ (σI) 2 τ σI + 1 + k r (21)
β(I) = k d τ (σI) 2 τ σI + 1 , (22)
with initial condition [START_REF] Han | Photosynthesis-irradiance response at physiological level: A mechanistic model[END_REF].

PROOF. It is straightforward to see that the equilibrium point  = 0 of the boundary layer model in ( 17) is asymptotically stable. By applying Thikhonov Theorem (see e.g. Theorem 11.1 in [START_REF] Khalil | Nonlinear Systems[END_REF]), we get that for any fixed T > 0 there exists ε 0 such that 0 < ε < ε 0 and system ( 11)-( 14) has an unique solution A t, ε , C t, ε on [0, εT ], and

C t, ε -C t = O(ε) (23) A t, ε - 1 -C t τ σ Ĩ t + 1 -A(t) = O(ε) (24)
hold uniformly for t ∈ [0, εT ], where C is the solution to equation ( 16) with initial condition [START_REF] Han | Photosynthesis-irradiance response at physiological level: A mechanistic model[END_REF]. Using the change of variable [START_REF] Nadarajah | A comparison of the continuous and discrete adjoint approach to automatic aerodynamic optimization[END_REF] we get that ( 18) and ( 19) hold for t ∈ [0, T ]. This concludes the proof of Lemma 1.

The net specific growth rate is given by the balance between photosyntesis and respiration µ(A, I) = kσIA-R where the constant R > 0 denotes the respiration rate and k > 0 is a coefficient. This gives the following dynamics for the biomass

X dX dt = µ(A, I)X . (25)
Now, as we are interested in the evolution of C we expressed the function µ as a function of this latter state.

Because A = 1-C στ I+1 , then the biomass dynamics writes

µ(A, I) = γ(I)C + ζ(I) = µ(C, I) , (26)
where

γ(I) = - kσI τ σI + 1 (27) ζ(I) = kσI τ σI + 1 -R . (28)
In the rest of the paper we assume that the light signal I is given by the function

I(t) = - 4 T 2 (I 0 -100) t 2 + 4 T (I 0 -100) t + 100 , (29)
where t ∈ [0, T ] and I 0 is the maximal light intensity of the signal. The parameters used are given in Table 1.

I0 1500 [µmol.m -2 .s -1 ] kr 4.8e -4 [s -1 ] k d 2.99e -4 [-] τ 1/0.1460 [s] σ 0.0029 [m 2 .(µmol) -1 ] k 3.6467e -4 [-] R 0.05 [d -1 ] Table 1
Parameter values used for the simulations.

PDE Modelling

The dynamics modelling of the RAB during the day can be represented by Partial Differential Equations (PDEs). It consists in decoupling the system in two: an enlightened part and a dark part. For each part we consider two PDEs: a PDE for the photoinhibition dynamics and a PDE for the biomass dynamics. Specifically, one has the following dynamics for the photoinhibition rate and the biomass, for the enlightened portion of the RAB,

∂ t C 1 + u(t)∂ x C 1 = -α (I(t)) C 1 + β (I(t)) (30) ∂ t X 1 + u(t)∂ x X 1 = γ (I(t)) X 1 C 1 + ζ (I(t)) X 1 , (31)
respectively and (t, x) ∈ R + ×(0, ), where is the length of the belt exposed to light (see Figure 1).The functions α, β, γ, and ζ are defined in ( 21), ( 22), ( 27), and (28) respectively. The light signal I(t) depends on the time because in outdoor cultivation photons come from the sun. In the above equations u denotes the velocity of the RAB. We suppose that the velocity can vary during the process operation and that it can be adapted to the state of the system together with the light intensity. For the portion of the RAB in darkness we have

∂ t C 2 + u(t)∂ x C 2 = -k r C 2 (32) ∂ t X 2 + u(t)∂ x X 2 = -RX 2 , (33)
where (t, x) ∈ R + × (0, (N -1) ). Here, N is the total length of the conveyor (1/N : ratio of belt which is enligthen, called the belt folding ratio). The space derivative corresponds to the transport part and the righthand sides of (30)-(33) correspond to the reaction terms. Since we are working in bounded space domains (0, ) and (0, (N -1) ) we need a boundary condition for each equation to get a well-posed problem. It is given by

C 1 (t, 0) = C 2 (t, (N -1) ) , C 2 (t, 0) = C 1 (t, ) (34) X 1 (t, 0) = X 2 (t, (N -1) ) , X 2 (t, 0) = X 1 (t, ). ( 35
)
The above boundary conditions may be read as "what leaves the enlightened part enters the dark part and vice versa". Finally, initial conditions shall be specified

C i (0, •) = C 0 i ∈ E , X i (0, •) = X 0 i ∈ L ∞ ((0, 1); R + ) , (36)
for i = 1, 2. We are interested to maximize the production of biomass between the intial time instant and t = T and we get the productivity by ponderating it by the production duration (productivity per unit of time). Therefore, we state the average areal productivity of our process as

P(u, N ; T, , I) = 1 T (N -1) 0 X 2 (T, x) -X 0 2 (x) dx + 0 X 1 (T, x) -X 0 1 (x) dx . ( 37
)
Proposition 2 The average areal productivity over the time interval [0, T ] can be rewritten as 31) and (33), the quantity P in (37) rewrites as

P(u, N ; T, , I) = - R T T 0 (N -1) 0 X 2 (s, x)dxds + 1 T T 0 0 (γ (I(s)) (X 1 C 1 ) (s, x) +ζ (I(s)) X 1 (s, x)) dxds . ( 38
) PROOF. Since T 0 ∂ t X i (t, x)dt = X i (T, x) -X i (0, x), using (
P = - 1 T T 0 u(t) (N -1) 0 ∂ x X 2 (t, x)dxdt + T 0 u(t) 0 ∂ x X 1 (t, x)dtdx + 1 T -R T 0 (N -1) 0 X 2 (t, x)dxdt + T 0 0 µ (C 1 (t, x), I(t)) X 1 (t, x)dxdt . ( 39
)
Using boundary conditions (35) and the expression of µ given in (26) one gets (38).

4 Gradient-Based Optimization of the RAB

Statement of the Optimization Problem

In this section, the objective is to take benefit of the model to optimize process performances and especially productivity as given by the equation in (38).

We use a practical method for an efficient identification of the optimal conditions. More precisely, we will find these conditions with an adjoint-based optimization algorithm. We aim at optimizing both the velocity u of the RAB and the shape parameter N which relies on the belt folding ratio. The mixed parameter optimization and optimal control problem for the given time horizon [0, T ], given light signal I, and a given , reads max

(u,N )∈U ×D P(u, N ; T, , I) , (40)
where D = [2, +∞) and U is a precompact set of functions, which will be specified later. The main difficulty in optimizing the productivity P is the presence of the parameter N in the bound of the second integral in (38).

To get ride of this issue we consider the functional

P(u, N ; T, , I) = - R(N -1) T T 0 1 0 X 2 (s, x) dxds + 1 T T 0 1 0 γ (I(s)) X 1 C 1 (s, x) +ζ (I(s)) X 1 (s, x) dxds , (41)
where X 1 , X 2 , C 1 , and C 2 satisfy the following system

∂ t C 1 + u(t) ∂ x C 1 = -α (I(t)) C 1 + β (I(t)) (42) ∂ t X 1 + u(t) ∂ x X 1 = γ (I(t)) X 1 C 1 + ζ (I(t)) X 1 (43) ∂ t C 2 + u(t) (N -1) ∂ x C 2 = -k r C 2 (44) ∂ t X 2 + u(t) (N -1) ∂ x X 2 = -RX 2 , (45)
in (t, x) ∈ (0, ∞) × (0, 1) subject to boundary conditions

C 1 (t, 0) = C 2 (t, 1) , C 2 (t, 0) = C 1 (t, 1) , (46)
X 1 (t, 0) = X 2 (t, 1) , X 2 (t, 0) = X 1 (t, 1) , (47)
and initial conditions

C i (0, •) = C 0 i , X i (0, •) = X 0 i , (48)
where

C 0 i (x) = C 0 i ( x) , X 0 i (x) = X 0 i ((N -1) x) . (49
We have the following lemma.

Lemma 3

The mixed parameter optimization and optimal control for the nonlinear infinite dimensional problems (40) and (50) are equivalent.

PROOF. As equation (49) leave to believe we will use a change of variables, namely

x1 = x 1 , x2 = x 2 (N -1) , (51)
C i (t, xi ) = C i (t, x i ) , X i (t, xi ) = X i (t, x i ) , (52)
for i = 1, 2. By computing the partial derivatives w.r.t. t and x of C i and X i , one has that system (30)-(33) becomes (42)-(45) through relationships (51)-( 52). Then, it can be verified that the productivity in (38) is given by (41) for the dimensionless system (42)-(48). It concludes the proof of Lemma 3.

While the optimization problem (40) needs some elaborate analysis, the optimization problem (50) is rather simple to solve in comparison. Let us mention that system (42)-( 48) is well-posed. The proof is detailed in [START_REF] Lamare | Gradientbased optimization of a rotating algal biofilm process[END_REF] and consists in using the characteristics method and noting that system (42)-(48) has a cascade structure: C 1 influences X 1 , C 2 influences X 2 . From now on and for the sake of simplicity we will drop the bar in the notation of C i and X i , i = 1, 2, and the dependence in the time variable whenever it is unnecessary.

Existence of an optimal control

The velocity is supposed to be time-varying. Thus, the optimization problem consists in finding an optimal point in an infinite-dimensional space. In our case this functional space is

U := C 1 (0, T ) |u min ≤ u(t) ≤ u max .
This control space ensures that, for each fixed N , there exists an optimal control for our problem. More precisely, we have the following theorem.

Theorem 4 Let N be fixed. There exists an optimal control u * in U maximizing the function J .

The proof of Theorem 4 is detailed in [START_REF] Lamare | Gradientbased optimization of a rotating algal biofilm process[END_REF] and consists in writing the solution of our system with a solution operator. This solution operator depends on the velocity u. Using the compactness of U and by the Arzelà-Ascoli Theorem any sequence (u β ) in U has a convergent subsequence (u β ). By considering the associated sequence of solution operator U β , it may be shown that the set of generators C is strongly compact in the sense of [START_REF] Slemrod | Existence of optimal controls for control systems governed by nonlinear partial differential equations[END_REF]. Then, by applying Theorem 3.1 and Theorem 5.1 of [START_REF] Slemrod | Existence of optimal controls for control systems governed by nonlinear partial differential equations[END_REF] we get that there exists an optimal control u * in U maximizing the function J .

The mixed optimization problem is generally numerically solved with a gradient descent algorithm. The derivation of the gradient of the cost function with respect to the velocity and the folding of the RAB cannot be derived as for a vector space function. The socalled "adjoint-method" may be used [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: Theoretical aspects[END_REF]. This method is based on an adjoint-system which has to be solved backward in time. Notice the adjoint variables play the role of Lagrange multiplier associated with the optimization constraints. This method is strongly linked to the calculus of variations theory. The interest reader is reported to the introductory monograph [START_REF] Liberzon | Calculus of Variations and Optimal Control Theory: A Concise Introduction[END_REF]. Finally, the continuous gradient is expressed as a function of the adjoint-state and nominal state solutions.

From now on we rewrite the velocity u as the sum of a positive constant term and an unknown function:

u(t) = u 0 + δ(t) , (53)
where u 0 ∈ R + and δ ∈ C 1 ((0, T ); R) such that u ∈ C 1 ((0, T ); R + ). The term u 0 turns out to be a key control which is more efficiently managed when explicited in (53) (see Remark 7 below). Besides, u 0 imposes the rotation direction.

Remark 5

The velocity is assumed to be continuously differentiable to ensure that the optimal control problem as a solution. But, in practice it would be continuously differentiable. Indeed, the velocity is induced by an engine and the inertia of the system leads to an equation of the form m u = F , where the right-hand side term F represents the balance between engine couple and friction terms, both terms being continuous.

Adjoint-Based Gradient Derivation

In order to consider a minimization problem we will consider the cost function in (41) times minus one i.e.

J (u 0 , δ, N ) = -P . ( 54
)
We have the following result.

Proposition 6 The gradient of the cost function, stated in (54), with respect to U = (u 0 , δ, N ) is given by

∇ u0 J (U ) = - 1 T T 0 1 0 p 1 ∂ x C 1 + q 1 ∂ x X 1 + p 2 ∂ x C 2 M + q 2 ∂ x X 2 M dxds (55) ∇ δ J (U ) = - 1 T 1 0 p 1 ∂ x C 1 + q 1 ∂ x X 1 + p 2 ∂ x C 2 M + q 2 ∂ x X 2 M dx (56)
∇ N J (U ) = u M 2 T T 0 1 0 (p 2 ∂ x C 2 + q 2 ∂ x X 2 ) dxds + R T T 0 1 0 X 2 dxds , (57)
where M = N -1 and (p 1 , p 2 , q 1 , q 2 ) satisfy the following system

-∂ t p 1 - u ∂ x p 1 = -α (I) p 1 + γ(I) (q 1 -1) X 1 (58) -∂ t q 1 - u ∂ x q 1 = (γ(I)C 1 + ζ (I)) (q 1 -1) (59)
-∂ t p 2 - u (N -1) ∂ x p 2 = -k r p 2 (60)
-∂ t q 2 - u (N -1) ∂ x q 2 = -Rq 2 + R(N -1) , (61)
with the boundary conditions

p 1 (t, 1) = p 2 (t, 0) N -1 , p 2 (t, 1) = (N -1)p 1 (t, 0) , (62) q 1 (t, 1) = q 2 (t, 0) N -1 , q 2 (t, 1) = (N -1)q 1 (t, 0) , ( 63
)
and the initial conditions

p 1 (T, •) = p 2 (T, •) = q 1 (T, •) = q 2 (T, •) = 0 . (64)
PROOF. For U = (u 0 , δ, N ) ∈ R + ×C 1 ((0, T ); R + )× R + , we denote by C i and X i , i = 1, 2, the corresponding solution to system (42)-(47). For a perturbation U = ũ0 , δ, N ∈ R + × C 1 ((0, T ); R + ) × R + we denote by C i and X i , i = 1, 2, the solution to system (42)-( 47) for the variables U + ε U . Denoting

C i = lim ε→0 C i (t, x) -C i (t, x) ε X i = lim ε→0 X i (t, x) -X i (t, x) ε J = lim ε→0 J U + ε U -J (U ) ε ,
we have

J = - 1 T T 0 1 0 γ (I) X 1 C 1 + γ (I) X 1 C 1 +ζ (I) X 1 dxds + R(N -1) T T 0 1 0 X 2 dxds + R N T T 0 1 0 X 2 dxds .
The variables C i and X i satisfy the so-called tangent model. Due to the space limitation, we only give the expression of the system satisfied by C 1 and C 2

∂ t C 1 + u 0 + δ ∂ x C 1 = - ũ0 + δ ∂ x C 1 -α (I) C 1 ∂ t C 2 + u 0 + δ (N -1) ∂ x C 2 = - ũ0 + δ (N -1) ∂ x XC 2 + N u (N -1) 2 ∂ x C 2 -k r C 2 ,
where ũ(t) = ũ0 + δ(t), with the associated boundary conditions

C 1 (t, 0) = C 2 (t, 1) , C 1 (t, 1) = C 2 (t, 0) , (65)
and the initial conditions

C 1 (0, •) = C 2 (0, •) = 0 .
Next we multiply the equations for C 1 , C 2 , X 1 , and X 2 by some functions p 1 , p 2 , q 1 , and q 2 respectively, representing the adjoint variables (see, for instance, [START_REF] Curtain | An Introduction to Infinite-Dimensional Linear Systems Theory[END_REF] page 601, for the precise meaning of the adjoint system of an infinite dimensional system), then we integrate by part over [0, T ]×(0, 1) and we obtain that p 1 , p 2 , q 1 , and q 2 satisfy system (58)-( 61) with the boundary conditions ( 62)-( 63) and the initial conditions (64) (for t = T since the system (58)-(61) has to be considered backward in time). Now, using that

lim ε→0 J U + ε U -J (U ) ε = ∇J (U ) • U ,
we get the formal expression for ∇J (U ) given by ( 55), (56), and (57). This concludes the proof of Proposition 6.

Remark 7 Let us remark that for t = T , the initial conditions (64) for the adjoint system yields

1 0 p 2 (T, x)∂ x C 2 (T, x) (N -1) + q 2 (T, x)∂ x X 2 (T, x) (N -1) dx + 1 0 (p 1 (T, x)∂ x C 1 (T, x) + q 1 (T, x)∂ x X 1 (T, x)) dx = 0 .
Hence, without the constant term u 0 in (53) we would not be able to control the final velocity of the rotating biofilm.

Besides, with constant initial conditions of the PDEs the same observation holds for the starting velocity.

Optimization Algorithm

In practice, we applied the following strategy to initialize and stop the gradient-based minimal search. In the numerical implementation of the descent algorithm we implement the gradient of J for the full-discretized system (discretization in time and space of the equations). There are two reasons guiding the full-discretization choice:

• It is a well-known fact that a space semi-discretization is not valid for a transport-reaction equation. Indeed, the derived equation get properties that the non-discretized equation does not have, such as the observability in any time while the non-discretized equation has only the property to be observable in a given minimum time T > 0. Therefore, the only possibility is a time semi-discretization. Hence, the resulting equation is a discrete-time system and not an ODE. Applying the maximum principle of Pontryagin (see [START_REF] Liberzon | Calculus of Variations and Optimal Control Theory: A Concise Introduction[END_REF]) can be difficult with the derived discrete time system. • But most importantly, at the end we have to discretize the equation for a practical implementation and it is well-known that it is more difficult to get stability results on the full-discretization than in the semidiscretization framework.

The strategy to derive the gradient in the continuous case is essentially the same as for the discrete case. In realistic operation u cannot reach any constant maximum value. Therefore, we impose an upper bound u max > 0 on the velocity. Besides, we also exclude the velocity u ≡ 0 which corresponds to the uncontrolled case. Hence, we impose a lower bound u min > 0. At each iteration, once the velocity is updated we proceed with a projection step of the computed velocity on the interval [u min , u max ], with the method proposed in [START_REF] Herzog | Algorithms for PDE-constrained optimization[END_REF].

The optimization program has been implemented in a C++ code. The Lax-Wendroff scheme (see e.g. [START_REF] Hirsch | Numerical Computation of Internal and External Flows: The Fundamentals of Computational Fluid Dynamics[END_REF]) is used to solve the transport-reaction PDEs (42)-(47). This scheme has been favored because it avoids a numerical diffusion in the solution and its accuracy is of second order in space and time. The number of points for the spatial discretization is 50 corresponding to dx = 0.02. The time step dt is governed by a classical Courant-Friedrichs-Lewy (CFL) condition ensuring the stability of the discrete scheme. This condition dictates a direct relation between the spatial mesh and the time mesh. By refining the spatial discretization, the time discretization will be thiner, giving a higher computation time.

5 Numerical Results

Process Exploitation Hypotheses

The set of parameters given in Table 1 are such that for an algae permanently exposed to a PFD larger than 168 µmol.m -2 .s -1 , the growth rate µ decreases when the value of irradiance I increases. This situation will be our reference case, and the gain in productivity with the RAB process will be computed in reference to it. The belt length exposed to light is fixed to = 1 m. Let us recall that the light signal considered in this paper will be the one given in (29). Therefore, N will give the total length of the belt in meters (we recall that it is given by N ).

The final time, to assess the productivity at the end of the day, is 28800 seconds (8 hours). The upper bound for the velocity is 0.3 m.s -1 , to reduce energy need for circulating the belt, and the lower bound to be 0.005 m.s -1 to avoid belt immobilization, and thus ensure nutrient supply to the microalgae. The tolerance of the gradient algorithm is fixed to 10 -4 .

Optimization of the RAB

The simulations were run on a 2012 commercial laptop with 4 Gb of RAM and a 2.5 GHz Intel Core i5 processor. First, we optimize the process with homogeneous initial conditions, namely,

X 0 1,2 = 10 [g/m] , (66)
C 0 1,2 = κ , (67)
where κ ∈ [0, 1]. The starting point for the optimization is (u, δ, N ) = (0.01, 0, 3 ). In other words, the optimal rate N opt provides an average irradiance received by the cells, close to the optimal one obtained at equilibrium in a fixed biofilm permanently exposed to light (N = 1). Thus, the RAB provides an optimal way to dilute light along time.

For every initial conditions tested, the optimal velocity computed is constantly equal to u max (0.3 m.s -1 ). We have carried out other simulations with a higher u max and we observe the same behavior. Besides, this observation is in accordance with the flashing effect already reported in the literature [START_REF] Hartmann | The effect of photosynthesis time scales on microalgae productivity[END_REF], demonstrating that the faster light varies, the higher the growth rate.

Remark 8 From an optimal control point of view the saturation of the velocity on its upper constraints looks like a bang-bang control. Unfortunately, there does not exist any result similar to the Pontryagin's Maximum Principle (PMP) for nonlinear systems in infinite dimension [START_REF] Egorov | Necessary conditions for optimal control in Banach spaces[END_REF].

To better understand the role of the velocity we run a simulation where the only parameter to optimize is N , the velocity u is fixed to 0.001 m.s -1 , the initial condition for the photoinhibition rate is 0.5. As shown by Figure A.2 the optimization process stops after 4 iterations, and N stays small compared to the case where the velocity is higher (cf Figure A.1). This result should be understood as follows: for a low velocity, having a large N means staying a long time in the dark and thus a large loss of biomass by respiration. Therefore, the velocity should be above a certain threshold. Once the velocity is large enough, the most impacting parameter is N which mainly drives the optimal solution. Finally, the adaptation of the velocity does not provide a significant gain, and it does therefore not seem mandatory.

Moreover, in order to minimize the energy input in the RAB process, we computed the (time-varying) lowest velocity for which the deviations with respect to the op-timal productivity is no more than 5%. The result is depicted on Figure 2. As we can see, the average value of this "energy saving velocity" is lower than the constant velocity computed from optimization process. It means that we may find some velocity controls for which the productivity remains close to the optimal productivity and which is less energy consuming. Adding a penalty term for the velocity in the cost function (41) would also result in decreasing the optimal velocity.

So far, we have supposed that the light signal I is time dependent and exogenous to the system, since the photons come from the sun. Now, we led some numerical experiments with a spatial dependent light signal:

I(x) = I 0 + A sin (ωx) . (68)
We have observed that the influence of the three parameters I 0 , A, and ω is of second order compared to N and u. For instance, for fixed N and u the RAB is less efficient with I given by (67) than with I ≡ I 0 .

Then, we have considered the situation of nonhomogeneous initial conditions for the biomass density and the photoinhibition rate:

X 0 1 (x) = 0.24 cos (2xπ) + 6.6 (69) X 0 2 (x) = 0.24 cos (2(x + 1)π) + 6.6 (70) C 0 1 (x) = 0.25 cos (2xπ) + 0.25 (71) C 0 2 (x) = 0.25 cos (2(x + 1)π) + 0.25 . (72)
Here, the starting point for the optimization is (u, δ, N ) = (0.1, 0, 3) . As for the non-homogeneous initial conditions case the optimal velocity computed is constantly equal to the upper bound u max . The evolution of the cost function and of the ratio N are depicted on Figure A.3, respectively. The ratio is smaller than in the homogeneous initial conditions case. This is due to the fact there is a proportion of the biofilm with a small density of biomass.

Finally, we have displayed the gain between the considered rotating microalgal biofilm process and a biofilm permanently exposed to the light (N = 1) on Figure 3 in function of different I 0 in (29). More precisely, we have computed the productivity for different initial photoinhibition rate conditions at each I 0 . Then, we have considered the average value of these productivity. As we can see the gain is always positive meaning that the investigated technology gives better productivity than a biofilm permanently exposed to light. Besides, the gain increases as the maximum light intensity becomes large. For instance, for I 0 = 2500 µmol.m -2 .s -1 the gain is greater than 100 %. Moreover, on the same figure, the average N for the different initial photoinhibition rate conditions at each I 0 is depicted. As expected, greater is the maximum light intensity greater is the optimized N . Indeed, a high maximum light intensity implies much more inhibition for photosystem PSII than at a low maximum light intensity. Thus, much more time in the dark is needed by the damaged photosystems to recover, meaning a larger N should be consider for the belt.

Conclusion and Perspectives

We have explored a new process to enhance microalgae productivity in the objective of reducing costs and environmental impact for the production of biofuel. We have demonstrated that the process acts as a light diluter in time, in order to reach an average light (per cycle) close to the optimal constant light for photosynthesis.

This study also highlights the power of an adjoint-based gradient method to identify the optimal working mode of a complex process described by PDE. This approach can be extended to a more general case, where algae are not fixed on a conveyor, but are advected in a classical raceway cultivation process. The full-discretization of the equations is suitable for the 3D shape optimization also because it becomes straighforward to account for the problem constraints (for instance, possitivity or mass conservation). Then, cost-functions depending on the energy consumption should be considered in the future. In order to treat this perspective a deepeer analysis should be led to derive a relevant criterion to penalize the energy demand in the biological production.

In conclusion, this class of problems has many industrial perspectives with a variety of applications for renewable resources. Theoretical applications developed for these systems are still in their infancy and we are convinced that solving this problem is an important step for the improvement of microalgae based processes.

Fig. 1 .

 1 Fig. 1. Schematic representation of the RAB. In this case the belt folding ratio is 1/3.

)

  And let us consider the new objective max (u,N )∈U ×D P(u, N ; T, , I) .

Fig. 2 .Fig. 3 .

 23 Fig.2. Velocity found when allowing a deviation of 5 % with respect to the optimal productivity for an initial condition C 0 i = 0.3.

Fig. A. 1 .

 1 Fig. A.1. Evolution of the cost function J in (54) (solid lines) and of the parameter N (dashed lines) with respect to the iteration for different homogeneous initial conditions.

Fig. A. 2 .Fig. A. 3 .

 23 Fig. A.2. Evolution of the cost function J in (54) (solid line) and of the parameter N (dashed line) with respect to the iteration for homogeneous initial conditions when optimizing only N for a velocity equal to 0.001 m.s -1 .

  ) . The convergence results obtained when starting from different constant initial conditions κ are displayed on Figure A.1. It turns out that less than 15 iterations are necessary to converge towards an optimal solution. And it took around one hour and forty minutes CPU time for a run.As expected the productivity is better with a lower initial photoinhibition rate, see Figures A.1. The optimization process converges towards a N close to 7.5 and this for each initial conditions of the PDEs tested, see Figures A.1. Therefore, the results of the optimization procedure seem to be independant of the initial photoinhibition rate. The average light on the RAB can be assessed by dividing the averaged impinging light (1030 µmol.m -2 .s -1 ) by this optimal N opt , yielding 137 µmol.m -2 .s -1 . The growth rate at equilibrium for this averaged light is approximatively 2.79 d -1 . It is worth remarking that this value deviates only from 1.06 % from the optimum growth rate at steady state (µ opt ≈ 2.82 d-1

Acknowledgments

The authors are grateful to the Inria Project Lab Algae In Silico and the Ademe project PHYTORECOLTE for their support.

A Appendix: Numerical Results of the Optimization Process

In this section, we collect figures obtained from the optimization procedure for the examples of Section 5.


Pierre Le Denmat

Thomas Gargot

Mohamed Chetouani

Dominique Archambault

David Cohen

Salvatore M Anzalone

The CoWriter robot: improving attention in a learning-by-teaching setup

Keywords: Social robot, Co-Presence, Virtual Agents, Writing disorders

In this paper we compare three learning-by-teaching scenarios in which a robot, a virtual agent or a voice, guide users and provide them feedback aimed at improving their handwriting abilities. The presented system is part of an effort focused on assessment and remediation of dysgraphia and dyspraxia. Results show the performances and the limits of the robot on inducing co-presence and attention.

Introduction

Developmental coordination disorder (DCD or dyspraxia) is a neurological disorder that impairs the acquisition and the execution of coordinated motor skills [START_REF] Association | Diagnostic and statistical manual of mental disorders (DSM-5 R )[END_REF]. It can be usually associated with learning disorders and in particular with disgraphya, a deficency in the ability of writing, affecting its execution in legibility and speed. Nearly 6% of children between 5 and 11 years in France is diagnosed with dysgraphia. Researchers highlight the fundamental importance of early intervention of handwriting skills, as soon as possible in children' educational path. With an onset in the early developmental period, the assessment of dysgraphia is performed through a standard test called Concise Evaluation Scale for Children's (BHK) [START_REF] Hamstra-Bletz | Concise evaluation scale for children's handwriting[END_REF]. This assessment can be difficult due to its high costs and subjectivity. However, recently [START_REF] Asselborn | Automated human-level diagnosis of dysgraphia using a consumer tablet[END_REF] has been shown how information and communication technology can help doctors in tackling such issues, by rapidly performing semi-automated, detailed diagnosis.

In this context, social robotics can help both the assessment and the remediation of children with dyspraxia and, more in general, with handwriting difficulties. A team from EPFL [START_REF] Hood | The cowriter project: Teaching a robot how to write[END_REF] recently proposed a learning-by-teaching scenario in which children teach handwriting to a robot through a tablet (Figure ??). The robot is able to fake its handwriting skills, proposing on the tablet purposefully deformed letters that take in account children performances and errors. Robot's handwriting skills can improve according to the quality of children' examples. The idea is to create an emphatic link between the child and its 'protégé', the robot, stimulating motivation and commitment. Children will act as mentors who help their protégé robot in the handwriting activities, getting themselves practicing without even noticing.

The main elements of this learning-by-teaching scenario are the robot and the tablet. This last one is used to collect handwriting samples from users, as well as to present handwriting feedback. It is possible to speculate about the importance of the presence of the robot and its ability on stimulating concentration and engagement. The experiment reported in this paper explores the role in this scenario of such agent, comparing users performances in three different conditions: handwriting sessions with the CoWriter robot; handwriting sessions with a virtual agent; handwriting session with the tablet only, guided by a voice. The hypothesis is that the social robot would be able to elicit an higher attention than the virtual agent or the vocal guide.

Materials and Methods

In this experiment we compare three CoWriter-based scenario experiments: in a first one, according to its original implementation [START_REF] Hood | The cowriter project: Teaching a robot how to write[END_REF], a user will interact with Nao robot, from Softbank Robotics (CR); in a second experiment the user will interact with a virtual agent, a simulated Nao (CA); in a third experiment, the user will be guided by a voice (CV).

A total of 12 adults (7 males, 5 females) in age between 12 and 31 years old (mean=23.9; std=2.75) participated to the experiment. Each participant inter-acted in the three conditions in a randomized way. For each condition, 5 words where chosen. For each word, in turn, the agent or the voice proposed its handwriting sample through the tablet; then the human user proposed his correction (Figure 2). The process repeated for each word, until the user positively judged the last handwriting sample.

At the end of the experiment, a questionnaire, the "Networked Minds Social Presence Inventory" was administered [START_REF] Biocca | The networked minds measure of social presence: Pilot test of the factor structure and concurrent validity[END_REF]. The purpose of this questionnaire is the measuring of the perceived social presence during an interaction, whether is face-to-face between two people, in telecommunication or with a non-human agent. While the original questionnaire is composed by six dimensions, this work focuses only in four facets, keeping apart the affective perception as not the focus of this experiment [START_REF] Harms | Internal consistency and reliability of the networked minds measure of social presence[END_REF]:

-Co-Presence: the user's awareness of the presence of the interaction partner; -Attentional allocation: the perceived attention received by the partner as well as the attention allocated towards the partner; -Perceived message understanding: the bidirectional communication understanding between the partners; -Perceived behavioral interdependence: the of the mutual behavioral connection between the partners.

Each dimension is composed by 6 questions and use a Likert scale 1-7 as response.

Fig. 2: The learning-by-teaching scenario implemented by the system.

Experimental Results

A one-way ANOVA test within subjects, comparing the three condition (CR, CA, CV) has been accomplished on the questionnaire data, revealing an effect on the Co-Presence dimension (F 2,22 = 8.69; p = 0.002). As in Figure 3 (left), post-hoc tests highlight a significant higher score in the robot condition compared to the voice condition (CR > CV , p=0.030) or the virtual agent condition (CR > CA, p = 0.039), but no difference is highlighted between virtual agent and voice (p = 1.000).

The ANOVA test revealed an effect also in the dimension of Attentional allocation (F 2,22 = 6.50; p = 0.006). As in Figure 3 (right), a significant higher score is found in the robot condition compared to the voice condition (CR > CV , p = 0.027), but still no difference between CV and CA (p = 0.595) and between CR and CA (0.085) is found.

The ANOVA test was not able to reveal any effect in the Perceived message understanding dimension (F 2,22 = 1.16; p = 0.33) as well as in the Perceived behavioral interdependence dimension (F 2,22 = 1.56; p = 0.23).

Conclusions and Future Works

Results on the Co-Presence dimension highlight that participants feel more the physical presence of the robot than its virtual agent. It is surprising, however, to note the absence of difference between the virtual agent and the voice. Results on the Attention allocation highlight how the robot is able to elicit more compliance than the vocal guidance. It should be noted in this case the absence of difference between the virtual agent and the robot. Such results can also highlight the possibility of the agent of capturing too much the attention of the user, acting as distraction of the task. The absence of difference in the two other dimensions, Perceived message understanding and Perceived behavioral interdependence, is not strange due to the particular task chose. In particular, the presence of a physical robot or a virtual agent seems does not impact on the comprehension of the exchange. At the same time, it is possible to hypothesize that the interpersonal exchange is too simple to be impacted by the presence of the artificial agents.

It should be noted how the small number of participant and their age do not permit a generalisation of the obtained results to children. Also, questionnaires information should be reinforced by behavioral measures (head movements, body movements, ...) that will be interpreted as the engagement state of the user during the shared task [START_REF] Anzalone | Evaluating the engagement with social robots[END_REF].

Fig. 1 :

 1 Fig. 1: A use case example of the CoWriter robot: the user teaches the robot how to write.

Fig. 3 :

 3 Fig. 3: The results of the questionnaires: co-presence, on the left; attentional allocation, on the right.