Thierry Colin
email: colin@math.u-bordeaux1.fr

François Cornelis

Julien Jouganous

Jean Palussière

Olivier Saut

Patient specific simulation of tumor growth, response to the treatment and relapse of a lung metastasis: a clinical case

Keywords: Clinical data assimilation, Medical imaging, Partial Differential Equations, Tumor growth modeling

In this paper, a parametrization strategy based on reduced order methods is presented for tumor growth PDE models. This is applied to a new simple spatial model for lung metastasis including angiogenesis. The goal is to help clinicians monitoring tumors and eventually predicting its evolution or response to a particular kind of treatment. To illustrate the whole approach, a clinical case including the natural history of the lesion, the response to a chemotherapy and the relapse before a radiofrequency ablation is presented.

Background

Introduction

The metastatic disease to the lung is frequently encountered in patients with cancer whatever the primary location and it has been associated with poor prognosis. The incidence of such disease in patients who have died of an extrathoracic malignancy is reported to be of 20% to 54% [START_REF] Hu | Histological subtypes of solitary pulmonary nodules of adenocarcinoma and their clinical relevance[END_REF]. Nevertheless, limited pulmonary metastatic disease can now be successfully treated not only for palliative reasons. By controlling the primary tumor and in the absence of widely disseminated disease in many organs, the resection or the ablation of pulmonary metastases may prolong survival, improve the quality of life and, in some cases, ensure cure [START_REF] Zisis | The management of the advanced colorectal cancer: management of the pulmonary metastases[END_REF]. During the last decade, a better management of the metastatic disease to the lung has been achieved by the evolution of imaging, medical oncology, and surgical techniques. There has been improvements in CT imaging quality and scan time [START_REF] Prokop | Lung cancer screening: The radiologist's perspective[END_REF], as well as advances in the field of nuclear medicine and MRI [START_REF] Ohno | New applications of magnetic resonance imaging for thoracic oncology[END_REF] which can give more precise information on the location and extent of the disease. In particular, there has been widespread use of PET/CT for evaluating patients with metastatic pulmonary disease, which can early detect metabolically active metastatic disease [START_REF] Lin | Prognostic value of preoperative metabolic tumor volumes on pet-ct in predicting disease-free survival of patients with stage i non-small cell lung cancer[END_REF]. The targeted oncology treatments can achieve improved responsiveness and the resection of the lesion turns to be possible now with minimally invasive surgical techniques as well as percutaneous thermal ablations (see [START_REF] Simon | Pulmonary radiofrequency ablation: long-term safety and efficacy in 153 patients[END_REF][START_REF] Von Meyenfeldt | The relationship between volume or surgeon specialty and outcome in the surgical treatment of lung cancer: a systematic review and meta-analysis[END_REF][START_REF] Von Meyenfeldt | Local treatment of pulmonary metastases: from open resection to minimally invasive approach? less morbidity, comparable local control[END_REF]). In order to continue this trend and to improve the knowledge of the pulmonary metastatic disease, some authors advocated the development of new tools able to explore the first steps of metastatic implantation [START_REF] Wood | The role of the tumor-microenvironment in lung cancer-metastasis and its relationship to potential therapeutic targets[END_REF]. The majority of pulmonary metastases derive from cancer cells that enter the lungs through the pulmonary arteries and disperse in alveolar capillaries. More rarely, metastatic disease is the consequence of lymphatic spreading or development directly in the bronchial tree (see [10]). Most of these cancer cells are able to adhere to the endothelium of the capillaries. However, the cell survival is ultimately determined by local and oncology factors. Although most of these malignant cells do not survive [START_REF] Zisis | The management of the advanced colorectal cancer: management of the pulmonary metastases[END_REF], those that survive penetrate the endothelium of the capillaries and install in the pulmonary parenchyma, where they grow. Therefore, a better understanding of this tumor growth could widely have direct clinical applications. Cancer growth modeling aims at describing, understanding and predicting the evolution of tumors using numerical models. This science is constantly evolving, from the first cellular automata adapted to the microscopic scale to Ordinary Differential Equations describing the global dynamics of the tumor or Partial Differential Equations taking the spatial distribution of the cancer cells into account. If modeling can be used to help biologists understand the complex mechanisms of tumor growth, in this work, the motivation is to develop a tool for clinical oncologists to improve patient monitoring and eventually to predict the response to a treatment of some metastases to the lung. As we deal with medical images, the spatial information is important and partial differential equations seem to be the most adapted method to model tumor growth in our context. A true difficulty is to include in a model the biological interactions, as complex as they are, responsible for the tumor evolution without neglecting that the model must be relevant for in vivo applications, and its parameters need to be recovered from clinical images. Among the diversity of PDE models existing on this topic, we chose, as a starting point, the one developed in [START_REF] Colin | System identification in tumor growth modeling using semi-empirical eigenfunctions[END_REF] for the promising results obtained. The diffusion term for nutrient is replaced by a micro model of angiogenesis (see [START_REF] Carmeliet | Angiogenesis in cancer and other diseases[END_REF]) which seems more relevant at the medical imaging scale. The goal is to simulate the evolution of a tumor from a given patient so we had to define a calibration method to personalize the model. This method must be a good trade-off between accuracy and computation time. For that purpose, the approach chosen is based on a reduced order method named Proper Orthogonal Decomposition (POD) (see [START_REF] Kunisch | Galerkin proper orthogonal decomposition methods for parabolic problems[END_REF] or [START_REF] Kunisch | Galerkin proper orthogonal decomposition methods for a general equation in fluid dynamics[END_REF]). The way we address this problem is similar to that of [START_REF] Colin | System identification in tumor growth modeling using semi-empirical eigenfunctions[END_REF]: we first design a PDE model of tumor growth dealing with the spatial distribution of cancer cells with respect to time T (t, x, y). This PDE model contains some parameters that are patient-specific. The problem is therefore to determine for which parameters the numerical solution fits the data. Once these parameters have been estimated (we call this step the calibration of the model), one can perform a prediction for a longer time scale. The outline of this paper is the following: in the following paragraph, we present a clinical case of a single metastasis including natural growth, response to a treatment and relapse. The Methods section is devoted to the presentation of the PDE model of tumor growth that we use and the data assimilation technics. Then the results are discussed in a dedicated part.

A typical clinical case

In this paper, we will focus on a clinical case of a patient with an isolated metastasis in the lung from a primary bladder tumor. The history of the patient is presented on Fig. 1.

Figures 1a, 1b and 1c are CT-scans of the natural history of the lesion, i.e. without treatment. The patient has been seen for a radiofrequency ablation on 2008/12/10. At this point, the lesion was too large and the oncologist started a chemotherapy. Figures 1d and1e are CT-scans during the treatment while Fig. 1f shows the relapse a month after the end of the chemotherapy. A radiofrequency ablation was performed in July 2009. The questions that we want to answer in this paper are the following ones:

1 Is it possible with Fig. 1a and Fig. 1b to estimate Fig. 1c? 2 Is it possible with Fig. 1a, 1b, 1c and Fig. 1d to estimate the efficacy of the treatment Fig. 1e and the amount of the relapse Fig. 1f?

Methods

A mathematical model for lung metastasis

The model we use in this work is derived from the one described in [START_REF] Colin | System identification in tumor growth modeling using semi-empirical eigenfunctions[END_REF]. We here consider only one kind of cancer cells. The tumor microenvironment, and in particular the amount of nutrients available, is essential to explain its evolution. Consequently, instead of directly modeling the nutrient density, we use a very simplified angiogenesis model to take into account the process by which the tumor escapes the avascular stage. The biological hypotheses of the model are the following:

• One single type of nutrient is considered. Its concentration controls cellular proliferation and death. • The amount of nutrients available is proportional to the quantity of blood vessels in the tissue. The nutrient diffusion is not taken into account as it is not relevant at the scale we consider. • The only cellular motion we take into account is the passive transport due to volume variations caused by mitosis or cellular death. • Cancer cells are continuously switching from the proliferative to the quiescent phenotype.

Cell behavior

The tumor cell density is denoted by T and satisfies the following equation:

∂T ∂t + ∇ • (vT ) = (γ p -γ d )T, ( 1
)
where v is the transport velocity. The left part of this transport equation satisfies the mass conservation principle and the source term drives this evolution (proliferation and death) of the cells population. The coefficients for proliferation and death by hypoxia (lack of dioxygen), γ p and γ d , are detailed in (2) and (3). Above a given threshold of nutrient supply M th , cancer cells tend to proliferate whereas below this threshold, they starve to death. The use of a hyperbolic tangent in both γ p and γ d expressions amounts to smoothing and regularizing the threshold functions, and K is a smoothing constant fixed. These functions are given by:

γ p (M ) = γ 0 1 + tanh(K(M -M th )) 2 , (2)
γ d (M ) = γ 1 1 -tanh(K(M -M th )) 2 . ( 3
)
From the tumor cells density T and the field γ p we can have an insight on the proliferative and quiescente cells location. Indeed the cells are proliferative were γ p = γ 0 and quiescente elsewhere, so proliferative celles are given by the expression P = γp γ0 T whereas quiescent cells are given by Q = (1 -γp γ0 )T . In addition to cancer cells, we also take into account the healthy tissue whose density is denoted by S. We assume here that this healthy tissue is transported at the same velocity as the cancer cells but is globally neither proliferating nor dying. More precisely, cellular birth and death counterbalance each other at the timescales we consider so we have, for S, a transport equation without source term (4):

∂S ∂t + ∇ • (vS) = 0. ( 4
)
As we consider that the tissue is saturated, the sum of the two densities T and S equals to 1 which gives the value of S (5) (see [START_REF] Bresch | Computational modeling of solid tumor growth: The avascular stage[END_REF]):

S = 1 -T. (5)
Summing ( 1) and ( 4) we obtain an equation on v (6):

∇ • v = (γ p -γ d )T. (6)
This is not sufficient to determine the velocity. However, to close the system of equations (see [START_REF] Ambrosi | On the closure of mass balance models for tumor growth[END_REF]), we consider that the velocity v is obtained using a Darcy law in Eq.( 7): v is derived from a pressure or potential π in the tissue.

v = -∇π. (7)
This formulation amounts to saying that tumor cells are pushed out if they are proliferating or pulled in if they are dying. Here we could also use a Stokes equation to describe the velocity (see [START_REF] Bresch | Computational modeling of solid tumor growth: The avascular stage[END_REF]) but it complicates the model without improving significantly its accuracy or biological relevance.

Angiogenesis

At the end of the avascular stage of its development, the tumor reaches such a size that its direct environment is not able to supply enough nutrients to allow the tumor to keep on growing. At this point, cancer cells emit chemical signals such as Vascular Endothelial Growth Factor (VEGF) which result in the emergence of a tumoral neo-vasculature (see [START_REF] Billy | A pharmacologically-based multi-scale mathematical model of angiogenesis, and its use in analyzing the efficacy of a new anti-cancer treatment strategy[END_REF][START_REF] Carmeliet | Angiogenesis in cancer and other diseases[END_REF][START_REF] Mantzaris | Mathematical modeling of tumor-induced angiogenesis[END_REF]). Research on this phenomenon is very active in particular to develop anti-angiogenic treatments. As most of the metastases visible on medical imaging are in the vascular stage, it was important to take the angiogenesis process into account. However, it is out of reach to model this extremely complex phenomenon in details. We therefore choose to use a very compact model in order to achieve a good balance between the complexity of the model and the feasibility of the data assimilation (see next section). It is described by the last two equations ( 8) and ( 9). The scalar variable ξ is the total amount of pro-angiogenic agent such as VEGF. It is secreted by quiescent cells given by the expression Ω (1 -γp γ0 )T dω, Ω being the computational domain, and evacuated by the organism.

∂ξ ∂t = α Ω 1 - γ p γ 0 T dω -λξ. (8)
As we assume that the quantity of nutrient is proportional to the density of blood vessels in the tissue, we mixed these two notions in the unique variable M we shall call "vasculature". It can be seen as the environment ability to meet the needs of the tumor in term of nutrients, and could be related with the carrying capacity introduced in the Hahnfeldt ODE model [START_REF] Hahnfeldt | Tumor development under angiogenic signaling: a dynamical theory of tumor growth, treatment, response and postvascular dormancy[END_REF]. The vasculature M is damaged by proliferating tumor cells while angiogenesis is promoted by quiescent cells; it therefore obeys

∂M ∂t = -ηT M + βξ(1 - γ p γ 0 )T. ( 9
)
The model takes into account some important mechanisms involved in tumor growth such as proliferation, death or angiogenesis. Moreover, it is simpler than the previous one [START_REF] Colin | System identification in tumor growth modeling using semi-empirical eigenfunctions[END_REF] as we just have one sort of cancer cells and there is no transport or diffusion of the vasculature. This trade-off is made to keep the model biologically relevant yet simple enough to be parametrized.

Taking therapeutical effects into account

The model architecture makes easy to include different types of treatment. Chemotherapy effects can be simulated by adding a death term on Eq.( 1) as -δT .

In this case Eq.( 1) becomes:

∂T ∂t + ∇ • (vT ) = (γ p -γ d )T -δT. (10)
To fullfil the saturation assumption Eq.( 6) is modified as follows:

∇ • v = (γ p -γ d -δ)T. (11)
We can also take into account an anti-angiogenic drug. In this case, we define a therapeutical indicator through a function denoted by f as:

f (t) =    1 without treatment, < 1 under treatment.
If we consider that it inhibits the production of pro-angiogenic agent, Eq.( 8) can be modified as follows:

∂ξ ∂t = αf (t) Ω (1 - γ p γ 0 )T dω -λξ.
If the anti-angiogenic drug also inhibits the growth factor membrane receptors on the endothelial cells, we can add a corresponding term to the vasculature equation Eq.( 9):

∂M ∂t = -ηT M + βξf (t)(1 - γ p γ 0 )T.
We postpone the modeling of this effect of anti-angiogenic drugs to a future work.

Data assimilation technique

The idea is to use the model that has been described in the previous section to obtain a forecast of the evolution of the lung tumor presented in the Background. This model contains numerous parameters that have to be recovered using the two images in order to perform the prediction. Some appear explicitly in the equations, such as the scalars α, β, γ 0 , γ 1 , η, λ and M th . Others are implicit and imposed as initial conditions, such as the scalar field M (t = 0) and the scalar ξ(t = 0), and in the case where the patient is under treatment at least one additional parameter has to be determined. The goal of this section is to build a calibration method fast and accurate enough to recover an adequate set of parameters describing the tumor evolution.

Simplifying assumptions

As mentioned before, we have 8 scalars and one spatial field to identify. These quantities cannot be estimated by in vitro or in vivo experiments. Furthermore, the parameterization problem is ill-posed. We can fix ξ(t = 0) arbitrarily without loss of generality since the variations of ξ(t = 0) can be taken into account with parameters α, β and λ. For convenience, concerning the order of magnitude in our numerical code, we take ξ(t = 0) = 0.1. For M , the situation is more complex and we have no universal solution to propose. We choose to take M (t = 0) = 0.8×T +S. This means that the available quantity of nutrients is lower inside the lesion than outside and that initially these quantities are constant. Note that this property is not satisfied for t > 0 because of equation ( 9). This is one of the strengths of the model: it is possible to obtain an heterogeneous distribution of nutrients within the tumor during the evolution that accounts for complex evolution and that makes the difference with scalar models dealing only with the volume of the tumor. Again the ratio 0.8 in the initialization of M is arbitrary and has to be related to M th and K.

Sensitivity analysis

In order to restrict the parameter space, we perform a sensitivity analysis on the model. It aims at determining if some parameters have more influence on the results than others. In this case, it means that the variation ranges of the less meaningful parameters must be modified or that we can simply fix them to a nominal value. For each parameter (p i ) i=1...7 = (α, β, γ 0 , γ 1 , λ, η, M th ), we choose a range of variation [a i ; b i ] given in Fig. 1. Denoting by T (t, x, y, p 1 , . . . , p 7 ) the solution given by the PDE system described previously, we compute numerically ∂ ∂pi Ω T (t f , x, y, p 1 , . . . , p 7 )dxdy for each i = 1 to 7, t f being the final time. We evaluate the mean value m i and the standard deviation S i of these quantities when the values of the parameters (p i ) i=1...7 are uniformly distributed on [a i ; b i ]. We use Morris method (see [START_REF] Morris | Factorial sampling plans for preliminary computational experiments[END_REF]) to perform this computation. The corresponding values of (m i , S i ) are plotted in Fig. 2.

The influence of the parameter can be established by the distance to the origin of the corresponding point (m i , S i ). In our case, we can see that the most influent parameter is γ 0 which is not surprising as it rules the exponential growth of the tumor. Conversely, the model is almost not sensitive to variations of α and β in the ranges we consider so a part of the inverse problem can be simplified by fixing them thus reducing the degrees of freedom.

Formulation of the inverse problem

Given a sequence of medical images or snapshots of the tumor, we aim at finding a suitable set of parameters able to reproduce the observed behavior. In other words, simulations with this set must fit as well as possible the clinical data. For this, our approach is to use a scalar objective function, which basically quantifies the distance between the observable data and the model simulation, and try to minimize it. There are different ways to measure this distance and as a criterion we chose a combination of the comparison of the mass and the L 2 distance as follows. The cost function is denoted by f obj :

f obj (p 1 , . . . , p 7 ) = w 1 × ns i=2 Ω (T data (t i ) -T model (t i , p 1 , . . . , p 7 )) 2 dω (12) + w 2 × ns i=2 Ω
T data (t i , x, y)dxdy -Ω T model (t i , x, y, p 1 , . . . , p 7 )dxdy , where (p 1 , ..., p k ) is the parameter set used in the model simulation, w 1 and w 2 are weights to balance the influence of each term, and n s is the number of snapshots we have. The first term in the objective function accounts for the shape while the second one only accounts for its volume. Of course, as initial data for our simulation, we use the "exact" value given by the data at time t = t 1 , therefore the sum starts at i = 2. As a consequence, we need at least two images of the metastasis at different times to personalize the model. The tumor cell density T data (t i ) is extracted from the i th snapshot (corresponding to time t i ). The quantity T model (t i , p 1 , ..., p k ) is the tumor cell density given by the simulation at time t i and for the parameters (p 1 , ..., p k ). To compute the L 2 distance, a data registration is necessary. We simply translate the images in order that their centers of mass match with the center of mass of the first image T 0 . The masses are computed by integrating the density on the domain Ω. From a computational point of view, whatever the minimization method is, one has to evaluate many times the objective function. It implies to simulate the model for lots of parameter sets which could be quite expensive. For instance, if we use a gradient algorithm to estimate seven parameters, at each iteration, the model is simulated eight times.

To make the calibration faster, we have developed a strategy based on a reduced order method called Proper Orthogonal Decomposition (POD) (see [START_REF] Sirovich | Low dimensional description of complicated phenomena[END_REF]).

Building a reduced order model to speed up computations POD resolution method for infinite dynamical systems consists in approaching partial differential equation systems with ordinary differential equations by decoupling the time and space variables (see [START_REF] Kunisch | Galerkin proper orthogonal decomposition methods for parabolic problems[END_REF] and [START_REF] Kunisch | Galerkin proper orthogonal decomposition methods for a general equation in fluid dynamics[END_REF]). The initial infinite dimension problem is replaced by a finite dimension (the smallest as possible) problem.

Let us first describe the POD used on the tumor cell density variable T . As we want to decouple space and time variables, we use the following representation for T (or any variable of interest):

T (X, t) = d i=1 a T i (t)Φ T i (X) + ǫ(X, t),
where a T i are scalar functions depending on time and Φ T i are spatial functions called modes and represent the geometry of the variable T . The dimension of the reduced problem is denoted by d. As it is an approximation, an error ǫ(t, X) is committed. The goal of POD is to provide us the best basis of spatial functions Φ T i to minimize the error. These functions Φ T i are extracted from a database of admissible behaviors of T . To generate this database, we sample the parameter space, run the direct model for each parameter set thus obtained and keep several snapshots (S T k ) k of the variable T . If the sample is correctly chosen, we have a representative set of geometrical configurations for the tumor cell density. Then we look for the functions Φ T i in the d-dimensional vectorial space generated by the snapshots from the database (we denote by E h = Span((S T k ) k ) this approximation space). These functions are taken as an orthonormal basis verifying: ∀h ∈ {1, ..., d}, ||T (X, t)-

h i=1 a T i (t)Φ T i (X)|| = min (Fi)i⊂E h ||T (X, t)- h i=1 a T i (t)F i (X)||.
This property is very interesting as it allows us to use a small number of modes without loosing too much accuracy.

The basis solution to this minimization problem is extracted from the autocorrelation matrix of the snapshots M T (M T i,j = (S T i |S T j ), (.-.) being a scalar product on Span((S T k ) k )) by the formula (13):

∀i ∈ {1, ..., d}, Φ T i (X) = 1 λ T i ns k=1 v T i [k]S T k , (13)
where λ T i denotes the i th eigenvalue of M T and v T i [k] the k th component of the i th eigenvector of M T . Once we have obtained these POD modes, we compute the projection of the PDE system on the approximation space E h using the eigenmodes Φ T i (X). We use the POD approach on both the tumor cell density T and the pressure field π which are the two fields driven by PDEs in our system. Thus we obtain the approximations

T (X, t) = n T i=1 a T i (t)Φ T i (X) and π(X, t) = nπ i=1 a π i (t)Φ π i (X)
, n T and n π being the numbers of modes we use to represent T and π.

Reduced order model on T:

We inject the POD expression of T in Eq.( 1), where we denote Γ = γ pγ d :

i ∂a T i ∂t Φ T i + i a T i ∇.(vΦ T i ) = i a T i ΓΦ T i .
Then this equation is projected on the orthogonal basis (Φ T k ) k , which gives:

∀j ∈ {1, ..., n T }, ∂a T j ∂t + i a T i (∇.(vΦ T i )|Φ T j ) = i a T i (ΓΦ T i |Φ T j ).
After simplification, we have a system of ODE on the POD coefficients (a T i ) i which is solved with an implicit Euler scheme:

∀j ∈ {1, ..., n T }, ∂a T j ∂t = i a T i (ΓΦ T i -∇.(vΦ T i )|Φ T j ). ( 14
)
That can be rewritten as:

dA T dt (t) = N T (t)A T (t),
where

N T i,j (t) = (Γ(t)Φ T j -∇.(v(t)Φ T j )|Φ T i )
. Therefore, we avoid the time interpolation phase that was essential in [START_REF] Colin | System identification in tumor growth modeling using semi-empirical eigenfunctions[END_REF] to approximate the time derivative. Consequently the hypothesis that two consecutive medical data are close in time is not necessary in our case and the accuracy of the scheme is improved.

Reduced order model on π: From Eq.( 7) and Eq.( 6) we obtain:

-∆π = ΓT.

Then we use the POD expression for π,

- j a π j ∆Φ π j = ΓT.
Finally we project this last equation on the orthogonal basis (

Φ π i ) i : ∀i ∈ {1, ...n π }, - j a π j (∆Φ π j |Φ π i ) = (ΓT |Φ π i ). ( 15
)
Here again, the system above can be written in the compact form:

N π A π (t) = -B π (t),
where N π i,j = (∆Φ π i |Φ π j ) and does not depend on the time and B π j (t) = (Γ(t)T (t)|Φ π i ).

Initialization: For the simulations, we also need to set initial values for the coefficients (a T i ) i and (a π i ) i ). These values are obtained by projecting the initial configurations of the fields associated on the POD basis: ∀i ∈ {1, ...n T }, a T i (t = 0) = (T (t = 0)|Φ T i ) and ∀j ∈ {1, ...n pi }, a π i (t = 0) = (π(t = 0)|Φ π i ).

Resolution of the inverse problem

The aim of this section is to sum up the algorithm we use in order to solve the inverse problem associated to the cost function Eq.( 12), that is find (p 1 , ..., p 7 ) = argmin(f obj (p 1 , ..., p 7 )). We proceed in two steps.

Step 1: We perform a Monte Carlo method in order to determine a first approximation of the parameters and to avoid a local minimizer. This consists in comparing simulations performed with parameters sets randomly chosen in an empirical parameters space. We keep the parameters set corresponding to the lower value of the cost function f obj (see Eq.12).

Step 2: Once we have this first approximation, we start a gradient descent method using the result of the first step as initialization. In this second step, the value of parameters α and β are fixed to those obtained in step 1, since the sensitivity analysis performed in paragraph 4.2 shows a low variability of the results with respect to their values.

Note that the use of the POD method implies that we deal with the evaluation of solutions of ODEs and not of PDEs during this process and this fact saves a lot of CPU time.

Results and discussion

Tumor growth The method previously described was applied to our clinical case (Fig. 1). The first scan (on 2008/06/07) represents the initial condition for the tumor cell density T . The other one (on 2008/09/22) is used in the objective function from Eq.( 12) as T data . The process to extract information from this type of medical images is quite simple. First, the tumor is delineated by a radiologist. The scan measures the permeability to X-rays which we assume is proportional to the density of cells. So we consider that the tumor cell density T is proportional to the normalized intensity of the corresponding voxel on the image. To further simplify the problem, we assume, as a first approximation, that tumor cells do not die from hypoxia so we fix γ 1 = 0. Let us remind that we also set the initial conditions for M and ξ. We still have six parameters to estimate. A parameter space is empirically defined and we generate a database of simulations on this space. In this case the database contains 374 simulations and for each simulation we save 5 snapshots for both T and π. From this collection of snapshots, we extract the POD modes and the calibration process is launched with the reduced model. For the simulation of the POD model we only use 36 modes on T and 10 on π which is here a good balance between accuracy and computation speed. We obtain the set of parameters given in Table 2.

These parameters are used to make a simulation on the direct model and try to predict the evolution of the tumor. A third CT scan on 2008/12/10 is available to compare our prediction and the real evolution.

The result of this simulation is compared with the observed growth of the metastasis on Fig. 3. The prediction from the model is satisfactory. To have a quantitative criterion of evaluation of the prediction, we can consider the mass of the tumor which is drawn on Fig. 4.

We can see that the horizontal error of prediction of the third scan is less than 10 days. Moreover, the relative error is below 5% of the total duration which is lower than the 10% segmentation error made by the clinician. If we look at the curves, we observe a small difference between the reduced model simulation and the direct model simulation. This is due to the POD hypothesis and the truncation of the POD basis. To quantify the accuracy in shape of the simulation, we calculated spatial indicators as the DICE and the volume concordance given by:

DICE = 100 × 2 * |T model T data | |T model | + |T data | and V C = 100 × (1 - |T model -T data | |T data |
). Moreover, the temporal prediction error is another significant indicator. Indeed, clinicians often need to know when the tumor will reach a critical size. It could be, for instance, the size over which the tumor cannot be treated using radio-frequencies. We denote by t i the time of the i th exam and t ′ i the time when the simulated tumor reaches the size of the real tumor at the i th exam. It is relevant to look at the delay between the simulation and the real case t i -t ′ i and the normalized delay 100× 3.

Chemotherapy

Given the fast growth of this metastasis, clinicians decided to treat the patient with a chemotherapy. The treatment starts just after the last scan on 2008/12/10 and ends on 2009/06/29. To monitor the efficiency of the treatment, three control scans were planned: two during the treatment (on 2009/03/21 and 2009/05/27) and a last one, one month after the end of the therapy (on 2009/07/27). This last one shows a relapse as the tumor starts growing again from the end of the treatment. It is interesting for us to calibrate the model on this decreasing phase to try to predict the response to the treatment and possibly the relapse. As the model has already been parametrized on the growth phase, we keep the same parameters from this calibration and use the cytotoxic drug modeling given in Eq.(10) and Eq. [START_REF] Colin | System identification in tumor growth modeling using semi-empirical eigenfunctions[END_REF]. This second inverse problem is easier as the only parameter to estimate is δ. This time again, we need two images. We use the last one before treatment (on 2008/12/10) as the initial condition and the first control (on 2009/03/21) to parametrize. We find δ = 1.05 and then the direct model is simulated up to after the date of the last scan (2009/07/27).

The evolution of the tumor mass during the treatment provides a good insight on the therapeutical efficacy. It is given on Fig. 5 and we can see that here again the model is predictive for this case and provides a good estimation of the response of the patient to chemotherapy.

The tumor shape during treatment obtained by simulation is compared to the clinical data on Fig. 6 and the shape indicators are given on Table 4.

The tumor shape is quite well reproduced by the model for the first control scan which we use to find the treatment parameter δ. However, on the last control scan made during the cure (on 2009/05/27), the tumor is quite hard to delineate due to various physiological phenomena that are not taken into account in the model (such as oedema or fibrosis). It results from the cytotoxic effects of the treatment on the tumor cells. Finally, the model thus calibrated gives a good prediction of the relapse after the end of the chemotherapy. For this last time point, shapes comparison provides a good result.

Another clinical test case

The whole calibration method described previously was used on another case of tumor growth. Here again, we use two scans at different time points to calibrate the model and a third image to quantify the accuracy of the prediction. The mass comparison is given on Fig. 7.

In this case, the growth is slower than in the previous case and the model is able to reproduce such a kind of dynamics. Indeed, the time error in prediction is about 5.6 days which, at the time scale we work and considering the tumor registration uncertainties, is a satisfactory result. The shape of the simulated tumor is also satisfying according to the spatial distribution of tumor cell densities on Fig. 8 and the indicators on Table 5.

Conclusion

To try to answer the questions raised by our clinical case a new approach was developed. First, a model was written that takes important mechanisms of tumor growth like hypoxia or angiogenesis into account. Yet the model is simple enough to be calibrated for a specific patient. We have described a new calibration method based on Proper Orthogonal Decomposition to fit the patient data and perform a prediction of the tumor evolution. Prediction which is confirmed qualitatively and quantitatively by the medical imaging. As the model contains various kinds of treatments we could, on the same clinical case, establish a precise quantitative prediction of the response to the treatment at the end of the protocol. As shown in Fig. 3, the evolution is really quick between 2008/09/22 and 2008/12/10 and is difficult to predict by the clinicians only using the scans. In that kind of circumstances, a mathematical approach like ours could have helped the oncologist in his diagnosis. On this particular case, the interest of the simulation is clear. For the time being we do not use any data on the real vasculature. With the constant progress of medical imaging one can imagine that such data will soon be available. This advancement will validate our modeling of the vasculature M . We mainly used here two dimensional data extracted from scans (that naturally give a 3D view of the tumor). The whole method was coded for three dimensional data, from the direct model simulation to the POD projection and the calibration algorithm. Currently, it is computationally too expensive to build a database in 3D so we cannot use the whole 3D process on the clinical case. However, by improving and making the database generation faster, it could be interesting to use a 3D method as the reduced model simulation calculation time is almost the same as in 2D. Moreover it would allow us to free ourselves from the choice of a particular slice. Table 1 Parameter space used for the sensitivity analysis.

ti-t ′ i ti-t0 , i = 1 , 2 .

 12 The values at the calibration time point (2008/09/22) and the final scan (2008/12/10) are given on Table

Figure 1

 1 FiguresTables

Figure 2

 2 Figure2Sensitivity analysis on the tumor growth model using the Morris method: effect of each parameter on the final tumor volume. The farther the point is from the origin, the more meaningful the corresponding parameter is for the model.

Figure 3

 3 Figure 3 Comparison between the prediction and the observation. Top raw shows the real medical image, below the simulation. The dates are on the left 2008/09/22 (calibration) and on the right 2008/12/10 (prediction).

Figure 4

 4 Figure 4 Mass evolution comparison between the data and the simulations for the tumor growth.

Figure 5

 5 Figure 5 Mass evolution comparison between the data and the simulations under treatment.

Figure 6

 6 Figure 6 Comparison between the prediction and the observation for the tumor decrease due to treatment and the relapse post treatment: above the real medical image, below the simulation. The dates are on the left 2009/03/21 (calibration) in the middle 2009/05/27 (prediction for the response to treatment) and on the right 2009/07/27 (prediction for the relapse).

Figure 7

 7 Figure 7 Second test case: mass evolution comparison between the data and the simulations.

Figure 8

 8 Figure 8 Comparison between the prediction and the observation for the second test case: above the real medical image, below the simulation. The dates are on the left 2010/03/11 (calibration) and on the right 2010/07/16 (prediction).

Table 2

 2 Parameter set obtained by the inverse problem.

	α	angiogenic agent synthesis factor	8.109
	β	angiogenic growth factor	7.241
	η	vasculature destruction factor	0.673
	γ 0	proliferation rate	1.116
	γ 1	death by hypoxia rate	0
	λ	angiogenic agent destruction rate	0.865
	M th	hypoxia threshold	1.045

Table 3

 3 Scalar indicators for the tumor growth of the first clinical case: DICE, Volume Concordance and delays.

		2008/09/22	2008/12/10
	DICE	90.96%	87.21%
	Volume Concordance	82.54%	77.76%
	Delay (days)	0	-6.7
	Normalized Delay	0%	-3.6%

Table 4

 4 Scalar indicators for the tumor under chemotherapy and rebound of the first clinical case: DICE, Volume Concordance and delays.

		2009/03/21	2009/05/27	2009/07/27
	DICE	92.26%	87.44%	84.79
	Volume Concordance	84.41%	74.56%	69.9%
	Delay (days)	0.5	-0.8	-6.4
	Normalized Delay	0.5%	-0.5%	-2.8%

Table 5

 5 Scalar indicators for the tumor growth of the second test case: DICE, Volume Concordance and delays.

		2010/03/11	2010/07/16
	DICE	85.41%	88.69%
	Volume Concordance	70.59%	76.45%
	Delay (days)	0	5.6
	Normalized Delay	0%	2.3%

Acknowledgements

This study has been carried out within the frame of the LABEX TRAIL, ANR-10-LABX-0057 with financial support from the French State, managed by the French National Research Agency (ANR) in the frame of the "Investments for the future" Programme IdEx (ANR-10-IDEX-03-02). Experiments presented in this paper were carried out using the PLAFRIM experimental testbed, being developed under the Inria PlaFRIM development action with support from LABRI and IMB and other entities: Conseil Régional d'Aquitaine, FeDER, Université de Bordeaux and CNRS (see https://plafrim.bordeaux.inria.fr/).

Competing interests

The authors declare that they have no competing interests.

Author details 1 Univ. Bordeaux, IMB, UMR 5251, F-33400, Talence, France CNRS, IMB, UMR 5251, F-33400, Talence, France Bordeaux INP, IMB, UMR 5251, F-33400,Talence, France INRIA Bordeaux-Sud-Ouest, F-33400,Talence, France. 2 Hôpital Pellegrin, CHU Bordeaux, Place Amélie Raba-Léon, 33000 Bordeaux, France. 3 Institut Bergonié, 229 cours de l'Argonne, 33000 Bordeaux, France.


Andrea Cattini

Devis Bellucci

Antonella Sola

Lech Pawłowski

Valeria Cannillo

Functional bioactive glass topcoats on hydroxyapatite coatings: Analysis of microstructure and in-vitro bioactivity

Keywords:

A bioactive glass topcoat was introduced to modify the surface morphology and in-vitro reactivity of hydroxyapatite (HA) coatings for biomedical applications. With this aim, a CaO-rich bioactive glass, termed BG_Ca (wt.%: 4.7 Na 2 O, 42.3 CaO, 6.1 P 2 O 5 , and 46.9 SiO 2 ), was selected due to its good bioactivity and low tendency to crystallize at high temperature. The standard HA coatings were sprayed through atmospheric plasma spray (APS) on steel substrates starting from commercial powders ("APS-HA" samples). The HA coatings, in turn, were subsequently coated with a thin layer of bioactive glass by suspension plasma spray (SPS), thus obtaining the duplex systems ("APS-HA/SPS-BG_Ca" samples). The samples with and without the BG_Ca layer were analysed by microstructural characterization and by in vitro tests in simulated body fluid (SBF). The analysis revealed an increased reactivity of the APS-HA/SPS-BG_Ca samples compared to the glass-free APS-HA coatings.

Introduction

Over the last decades, biomedical sciences have acquired a growing importance. In particular, the progressive ageing of population and diffusion of morbid obesity, with related musculoskeletal disorders, are driving the development of innovative materials for bone implants [1]. Though the production of biomedical devices is nowadays a mature and well-established practice, the research is constantly evolving with the aim of improving the reliability and effectiveness of the available prostheses, thus increasing their life expectancy in service conditions [2].

Hydroxyapatite, HA (Ca 10 (PO 4 ) 6 (OH) 2 ), is a calcium phosphate ceramic whose composition and structure are very similar to those of the mineral component of bone. It is therefore osteoconductive, with an excellent biocompatibility which allows bone cells to grow on its surface [3]. For this reason, HA is currently applied in orthopaedics and dentistry with very encouraging results [4]. Unfortunately HA is brittle and hence it is not suitable for structural applications. A way to overcome this problem is to use HA as a coating on a metal implant, which would provide the required mechanical strength. Nowadays, the plasma spraying is the standard method to produce HA coatings [5,6]. In fact, this method makes it possible to achieve a high deposition rate, which results in high coating production rate; moreover plasma sprayed HA coatings, in turn, may reach a suitable adhesion to the substrate [7]. Although the plasma spray process has been thoroughly studied and optimized in recent years, some isolated cases of implant failure are still reported and revision surgeries are sometimes necessary [8].

Moreover the host tissue-HA bonding reactions begin after several hours after implantation, and therefore HA coatings are not the optimum solution if a fast fixation is required [4,9]. Consequently, numerous studies address the development of new bioactive coatings [10][11][12]. Similarly, much research has been performed on both alternative materials and innovative deposition techniques [13][14][15][16]. An example of an emerging family of biomedical materials is given by bioactive glasses. They are particularly interesting, because they possess the highest degree of bioactivity controlled by their chemical composition [17]. In particular, 45S5Bioglass® (in wt.%: 24.5 Na 2 O, 24.5 CaO, 6 P 2 O 5 , 45 SiO 2 ) is the most representative glass of the so-called "Bioglass" family, first proposed by Prof. Hench at the end of the 1960s [18]. Due to their particular composition, bioglasses possess a strong bone-bonding ability; specific formulations are even able to connect to soft tissues [19]. Unfortunately, conventional bioglasses are likely to devitrify at high temperature and the development of crystalline phases is known to retard the original bioactivity [20]. Some new experimental compositions are about to emerge aiming at limiting the tendency to crystallize i.e. preserving the amorphous structure of the glass during hightemperature processes, such as sintering and thermal spray coating deposition [21,22]. A new strategy to obtain improved bioactive coatings may lead to the improvement of the deposition process. Indeed, even if plasma spraying is the most widely used technique to produce biomedical coatings, many innovations have been recently proposed, such as suspension plasma spraying (SPS) [23]. In this process, the feedstock is a liquid suspension of fine, nanometric or micrometric solids instead of the usual dry powder having micrometric size. This is a key point, since the dispersion of a solid powder in a liquid medium (water or organic solvent) makes it possible to process very fine particles or hygroscopic materials (such as bioactive glasses), which otherwise could not be employed in conventional thermal spraying because of their poor flowability. As a result of spraying process, a controlled microstructure, with very small (near-nanometric) intralamellar crystalline grains, can be obtained [24].

In the present study, a bioactive glass topcoat was applied by the innovative SPS technique on a conventional HA plasma sprayed coating. The BG_Ca glass (having composition, in wt.%: 4.7 Na 2 O, 42.3 CaO, 6.1 P 2 O 5 , and 46.9 SiO 2 ) was selected since it combines high bioactivity with a relatively low tendency to crystallize at high temperature [25,26].

The aim of such functional duplex system was to merge the HA's biocompatibility and long term stability, which are useful to protect the metal substrate from the physiological environment, with the high bioactivity of the SPS bioactive glass top-coat.

In order to assess the effectiveness of the bioactive glass topcoat to accelerate the mineralization process in a physiological environment, the samples with and without a BG_Ca layer were tested in vitro, by soaking them in a simulated body fluid (SBF) [27].

Materials and methods

Feedstock materials

The hydroxyapatite (HA) powder was a commercial product manufactured by Tomita (Japan). The powder had a mean diameter of d 50 = 120 μm, but it consisted of spray dried sub-micrometric particles having internal porosity of about 12% [28].

As described in a previous contribution [28], the bioactive glass (BG_Ca) [25,26] was produced by the conventional melt-quenching method. Appropriate raw materials were mixed and then molten in a platinum crucible at 1450 °C. The molten glass was poured into water. The obtained frit was dried overnight at 110 °C and dry milled in an agate mortar and sieved at 63 μm. The powder was further milled in an attrition system by using 0.8 mm zirconia balls and adding the dispersant Beycostat C213 (2 wt.% of the dry powder). Since bioactive glass powders are water-reactive [29], ethanol was selected as the suspension liquid. The obtained BG_Ca powder had a monomodal distribution with a mean diameter of d 50 = 4.7 μm (determined with a granulometer Partica LA-950V2, Horiba). The feedstock suspension was formulated by dispersing 20 wt.% of solid phase in 80 wt.% of ethanol.

Coatings' deposition

Plasmaspray runs were performed using a SG-100 torch (Praxair, S.T., Indianapolis, IN, USA). Substrates were 316L stainless steel disks (diameter: 25 mm; thickness: 8 mm), preliminary sandblasted to a roughness R a = 4.3 μm (determined by means of Tester T500, Hommel) and ultrasonically cleaned with ethanol and acetone. 316L steel was chosen because it has proved to be acceptable for orthopaedic implants over decades of use and currently it is widely diffused for fracture repair devices [START_REF] Pilliar | Biomedical Materials[END_REF].

During the APS deposition, the HA powder was supplied by a laboratory powder feeder including an Ar-pressurized chamber and a vibrating system.

To deposit the BG_Ca topcoat, a pneumatic system was used to feed the suspension (Fig. 1a); to avoid agglomeration and sedimentation, the suspension was continuously agitated by a magnetic stirrer. Both feedstock materials were (separately) supplied into the plasma by internal injection. The injection orifice was radial inside the anode and mounted at 16 mm from torch's exit. A schematic representation of the SPS deposition process is shown in Fig. 1b.

Several deposition attempts were carried out using tentative parameters taken from the literature [28,[START_REF] Dyshlovenko | [END_REF]. Then, based on the preliminary trials, the final spray parameters were defined, with the aim of optimizing the density and the thickness of the APS-HA layer and, at the same time, promoting the bioactivity of the glass topcoat, which should remain amorphous and properly porous. The deposition parameters set for the two powders are specified in Table 1.

During the spraying sessions, the temperature of the coatings' surface was measured using an optical IN 5 Plus pyrometer (Impac).

In vitro tests

In order to investigate the effect of the glass topcoat on the reactivity in physiological fluids, both the APS-HA samples and the APS-HA/ SPS-BG_Ca samples were soaked in a simulated body fluid (SBF) according to the standard procedure proposed by Kokubo [27]. The SBF was refreshed every 2 days to mimic the dynamic conditions of the human body. The samples were removed from the SBF after time periods of 1, 3, 7, and 14 days. The samples were washed with distilled water and left to dry at room temperature before further characterization.

Characterization

The samples were analysed before and after soaking in SBF by different techniques.

The surface of the coatings was investigated by X-ray diffraction, XRD, with an X'pert PRO diffractometer (PANalytical, Almelo, The Netherlands) equipped with the X'Celerator detector, and employing Cu-Kα radiation. Diffraction patterns were collected in the 10°-70°2θ range (step size: 0.017°, step time: 71.12 s).

An accurate microstructural analysis was carried out on the coatings' surfaces and polished cross sections by means of an environmental scanning electron microscope, ESEM (ESEM Quanta 200-FEI Company, Eindhoven, The Netherlands) coupled with an X-EDS microanalysis system (Inca, Oxford Instruments, U.K.). The microscope was operated both in high-vacuum and low-vacuum modes (pressure: 60 Pa).

The coatings' surface was analysed by micro-Raman spectroscopy (Jobin-Yvon Raman Microscope Spectrometer, HORIBA Jobin Yvon S.A.S., Villeneuve d'Ascq, France). The laser source emits 20 mW at 632.8 nm; 100× objective and no filters were employed.

Results and discussion

During the deposition of the APS-HA layer, the surface temperature did not exceed 250 °C and, at the end of the process, the temperature was about 200 °C (Fig. 2). During deposition of the SPS-BG_Ca topcoat, the temperature was slightly higher, with an average value of about 250 °C and local peaks of 330 °C. However such values are relatively low with respect to those usually achieved in SPS deposition with ethanol [28]. In fact, the final parameters (listed in Table 1) were purposely defined to avoid an excessive heat flux to the coating in order to prevent the glass from crystallizing. This also limited the sintering process and therefore it resulted in a porous glass topcoat, which is expected to enhance its bioactivity [28]. To reduce the heat flux to the sample, the spraying distance was set to 70 mm and the glass topcoat was sprayed in four discontinued spray runs, with two torch passes for each run. Moreover, in order to control the surface temperature, the samples were cooled down to the temperature below 70 °C before starting each runs (Fig. 2b). As shown in Fig. 3, the XRD confirmed that the HA coating basically preserved its phase composition after the deposition. In fact, the spectrum of the APS-HA coating is similar to that of the original HA feedstock powder. However some secondary phases could be also detected, which resulted from a partial decomposition of the HA. In particular, the peaks in the 29°-31°2θ range could be associated to the formation of tricalcium phosphate, TCP, and tetra-calcium phosphate, TTCP. The peak at 2θ = 37.5°corresponded to CaO. Since the relative amounts of the decomposition phases in HA coatings may be linked to the relative intensity of the corresponding peaks in XRD patterns [32], it is likely that the thermal decomposition of the APS-HA coatings was very limited, since the peaks associated to the secondary phases were relatively weak. The limited decomposition of HA could have been favoured by the relatively low temperature of particles in-flight in the plasma jet, resulting from the low enthalpy of the plasma gas. Indeed, it is known from the literature that HA is prone to decompose [33], since it loses its OH -groups (de-hydroxylation) at high temperature, with a twostep process resulting first in the formation of oxyhydroxyapatite and then of oxyapatite [34]. If the temperature does not exceed a critical point, these new phases can return into HA phase in the presence of water [35]. Otherwise, if the critical temperature is surpassed, the HA crystal structure changes and suffers a complete and irreversible dehydroxylation [36]. If the temperature increases further, HA starts to decompose to other calcium phosphates, such as TCP and TTCP. These phases, in turn, can degrade into calcium oxide and phosphorus pentoxide [34,36]. Even if the thermal decomposition is affected by several circumstances, such as the environment [37] and the initial purity and degree of crystallinity of HA [38,39], the critical point beyond which the decomposition is irreversible is about 1350 °C [36].

Fig. 3 also shows that, in the diffraction spectrum of the APS-HA/ SPS-BG_Ca sample, a broad band between 25°and 35°(2θ) was present, as often observed in silicate glasses. In addition, some new weak peaks could be detected, which may be associated with a calcium silicate, Ca 2 SiO 4 (larnite). It is interesting to note that in the relatively Na 2 Opoor BG_Ca glass a calcium silicate, instead of a sodium-calcium silicate developed, which is the typical outcome of the crystallization of 45S5 Bioglass® and other Na 2 O-rich bioglasses [26]. Moreover, well-defined peaks could be attributed to the HA layer underneath the glass topcoat. The partial transparency of the topcoat to the X-rays resulted from the intrinsic properties of silicate glasses, as previously observed for other silicate-based bioactive glass coatings [40], and from the structure of the topcoat itself, which was highly porous and relatively thin, as confirmed by the ESEM inspection described below.

The microstructure of the APS-HA samples was compact, with relatively low porosity, as shown in Fig. 4a. Some dark areas could be observed in the cross section. Since the EDS analysis (Fig. 5) confirmed that they were composed of HA, with atomic proportions analogous to those of the surrounding light-grey material, it is possible that such dark areas were the cross sections of unmelted HA particles incorporated into the coating. In fact, depending on the deposition conditions, partially molten or completely unmelted particles may remain in HA coatings produced by thermal spray techniques [41]. The surface of the APS-HA coatings (Fig. 4b) displayed both large splats, with diameters measuring some tens of microns, and areas with finer grains.

The inspection of the cross section of the APS-HA/SPS-BG_Ca samples confirmed that the glass topcoat was relatively thin, with typical thickness of about 20 μm (Fig. 6a). Most of all, the glass layer was characterised by a very fine microstructure, which is a distinctive feature of SPS coatings [23]. Moreover the BG_Ca topcoat was more porous than the atmospheric plasma sprayed HA layer. The rich porosity that characterises the glass topcoat is expected to enhance the implant osteointegration in vivo, since it confers a "trabecular-like" structure to the coating [42]. With regard to the surface morphology, the glass topcoat was extremely rough, as shown in Fig. 6b. A high roughness may be extremely advantageous, since it increases the specific surface area available for chemical reactions and hence favours the cell attachment [43]. Moreover, the surface roughness and morphology have significant impact on protein adsorption [44]. To conclude, very fine needle-like structures could be observed on the surface of the APS-HA/SPS-BG_Ca coatings. It is possible that the structure results from the crystals of the Ca 2 SiO 4 phase identified by the XRD.

After the immersion in SBF, the surface of both types of coatings (with and without the glass topcoat) was covered by dome-like structures which are typical of hydroxy-carbonate-apatite (HCA) deposits grown in vitro [45,46]. However, a direct comparison between the two types of samples (Fig. 7) suggests that the APS-HA/SPS-BG_Ca system reacted more rapidly than the APS-HA one. Indeed, the surface of the glass topcoat quickly assumed a more globular morphology. Moreover it is interesting to note that the fine crystalline structures originally present on the APS-HA/SPS-BG_Ca sample completely disappeared already after 1 day in SBF. Fig. 8 propose some representative images of the cross section of the samples soaked in SBF. After one day, a few micron-thick reaction layer grew on the APS-HA coating (Fig. 8a); as shown by the X-EDS analysis (Table 2), such layer had a Ca/P ratio of 1.49, a value observed already in calcium-deficient amorphous calcium phosphate (ACP) [47].

Instead, if the APS-HA/SPS-BG_Ca system is considered (Fig. 8b), after one day in SBF the glass topcoat underwent a through-thickness reaction and only silicon and oxygen could be detected by the X-EDS analysis (Table 2), thus suggesting a silica gel layer developed on the surface. The very high reactivity of the topcoat was due to the concomitant effect of the intrinsic bioactivity of the BG_Ca glass [26] and of the porous microstructure of the layer. In fact the rich porosity of the topcoat favoured the penetration of SBF in depth and greatly increased the specific surface area available for the chemical reactions. Moreover it was observed that, already after one day in SBF, the silica gel layer, in turn, was covered by a very thin deposit of calcium phosphate. The progressive conversion of the original glass to silica gel and the subsequent growth of the hydroxy-carbonate apatite (HCA) surface layer detected on the APS-HA/SPS-BG_Ca system are similar to the mineralization reactions occurring on conventional bioactive glasses in bulk form [48]. The growth of the calcium phosphate layer on the APS-HA/SPS-BG_Ca sample continued for increasing immersion time and, after one week in SBF, the original glass topcoat has disappeared, since it was completely converted to HCA (Fig. 8c) with a Ca/P ratio close to that of stoichiometric HA, i.e. Ca/P ~1.67 [3] (Table 2). It is interesting to note that the HCA precipitates were perfectly integrated into the coating, since the interface was not easily detectable.

As a general consideration, according to the SEM outcomes, the reaction kinetics proceeded faster on the APS-HA/SPS-BG_Ca system than on the pure APS-HA coating. In fact, the HCA layer grown on the APS-HA sample after 14 days in SBF was thinner than that grown on the APS-HA/SPS-BG_Ca system after 7 days in SBF, as apparent from the comparison of the corresponding SEM images in Fig. 8d andc.

The XRD pattern of the APS-HA sample reported in Fig. 9 (bold traces) did not show sensible variations as a result of the immersion in SBF. However, starting from the first day of immersion, the peaks associated to the secondary phases (TCP, TTCP and CaO) gradually disappeared. More precisely, the peak associated to CaO vanished already after one day of immersion, whereas the peaks of TTCP and TCP were not visible after seven days. On the other hand, the peaks of HA did not vary with increasing immersion time.

In contrast to this, the XRD patterns of the APS-HA/SPS-BG_Ca samples (Fig. 9, thin traces) strongly changed after immersion in SBF. Starting from the first day, the main peaks of HCA could be identified, although they were very broad. The large width of the peaks was indicative of the poor crystallinity of HCA, an effect that is often observed for in vitro grown HCA [49]. At the same time, the peaks due to Ca 2 SiO 4 disappeared, supporting the hypothesis that such crystallizations were mainly constituted by superficial fine crystals, readily dissolved by the SBF. The peak at 2θ = 29.4°, observed occasionally in the XRD spectra of the soaked samples, is due to CaCO 3 , a mineral that may co-precipitate in vitro with HCA [27,50]. Further information on the reactions occurring in SBF was obtained by means of Raman micro-spectroscopy. Fig. 10 compares the spectra of the APS-HA samples (bold traces) and those of the APS-HA/ SPS-BG_Ca ones (thin traces). In particular, if the pattern of the as-sprayed APS-HA coating is considered, the peak associated to the symmetrical stretching of the PO 4 group, at about 960 cm -1 , can be separated in two different contributions, namely a broad and weak peak, and a sharp and intense one. The latter can be associated to HA, while the former is thought to be due to the decomposition phases: ACP, TTCP and TCP [51]. After immersion in SBF (Fig. 10, bold traces), the Raman peaks associated with HA became broader and less defined, but they did not shift to different wavenumbers. The gradual change in the peak shape was likely due to the precipitation of poorly crystallized HCA on the original surface.

The Raman spectra acquired on the surface of the APS-HA/ SPS-BG_Ca coatings showed a marked change after the immersion in SBF (Fig. 10, thin traces). In the spectrum of the as-sprayed coating, three bands related to the BG_Ca glass were observed: a first band at about 635 cm -1 and a second one in the range from 800 to 1150 cm -1 , which are both associated to the Si\O\Si bonds in silica, and a third band at 950 cm -1 , which is due to the symmetrical stretching of the PO 4 group [52]. Then, after soaking in SBF, the new peaks of the precipitated HCA progressively replaced the signals of the original glass and they steadily became sharper with increasing immersion time. From this point of view, the evolution of the Raman spectra of the APS-HA/ SPS-BG_Ca samples in SBF corresponded to that of the XRD patterns, since in both cases the growth of HCA could be clearly detected.

Moreover, with regard to the Raman spectra, it is worth noting that, after immersion in SBF, a new peak appeared at approximately 1070 cm -1 , which is assigned to the carbonate group, thus confirming that the precipitated HA was indeed carbonated [52]. The comparison between the spectra of the APS-HA coatings and the APS-HA/ SPS-BG_Ca (Fig. 10) shows that the peak at 1070 cm -1 is stronger for the latter. Since the peak is strictly connected to the carbonate group of the HCA grown in SBF, this discrepancy between the Raman spectra further confirms that the APS-HA/SPS-BG_Ca samples reacted in SBF more quickly than the APS-HA ones.

Currently the discussion is still open on the real significance of in vitro tests. Indeed, the SBF (like other media, such as Hanks' balanced salt solution HBSS, Na-lactate and lactic acid-buffered physiological solution [53]) is simple salt solutions and, since it does not contain cells, it is not able to reproduce the real complexity of a biological environment [54]. However it is generally recognized that SBF tests may be useful to evaluate the ability to promote the growth of a HCA surface layer, which is similar to the mineral component of human bone [27]. So, even if the apatite forming ability does not automatically imply the bone-bonding ability [55], the development of a HCA layer in SBF may be still considered a favourable indicator of the tendency to support the mineralization process [27]. The results of the SBF tests therefore suggest that the addition of a bioactive glass topcoat may be effective to improve the mineralizing ability of a standard plasma sprayed HA coating, especially if the glass topcoat possesses a fine microstructure, with a controlled porosity.

Conclusions

In this study standard HA coatings with and without a functional bioactive topcoat were produced and tested. In particular, the HA coatings were sprayed through atmospheric plasma spray on stainless steel substrates, while the glass topcoat was produced via suspension plasma spray.

Table 2

Coatings' composition evaluated by X-EDS analysis performed on the cross-section of SBF soaked samples (Fig. 8). The bioactive glass topcoat showed high porosity and rough surface, which constitute an optimum microstructure to promote bioactivity. The coatings, with and without the glass topcoat, were tested in vitro by immersion in a simulated body fluid (SBF). It is interesting to observe that the samples with the glass topcoat showed a greater reactivity with respect to pure HA coatings. In fact, both types of coatings were covered by a calcium-phosphate layer after soaking in SBF, however the bioactive glass topcoat was completely converted into hydroxy-carbonateapatite in less than one week. This phenomenon was confirmed by both ESEM investigation (with X-EDS microanalysis), and by X-ray diffraction. Furthermore, the Raman analysis suggested that the HCA developed on the bioactive glass topcoat exhibited a greater carbonation level with respect to the layer that developed on pure HA.

Point

Nevertheless, this is only a preliminary trial, since cell tests are necessary to determine the biological effectiveness of the topcoat. In this way it could be possible to verify whether the high porosity of the topcoat, associated to the intrinsic bioactivity of the glass, could be able to increase the absorption of metabolites and the cellular activity.

Fig. 1 .

 1 Fig. 1. Coatings' deposition process: (a) SPS set-up, (b) production scheme of the coatings.

Fig. 2 .

 2 Fig. 2. Substrate temperature measured during spray sessions: (a) atmospheric plasma spray deposition of HA, (b) suspension plasma spray deposition of BG_Ca.

Fig. 3 .

 3 Fig. 3. Comparison between the X-ray diffraction pattern of the commercial HA feedstock powder and those of the coatings. Symbols: *, larnite (Ca 2 SiO 4 , 01-070-0388); ▲, TTCP (00-025-1137); ▼, TCP (00-009-0169); ■, CaO (01-075-0264); the peaks without symbol are due to HA (01-089-6437). Fig. 4. APS-HA coatings: (a) cross section, (b) surface.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Details of the APS-HA cross section and corresponding EDS microanalyses.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Surface morphology after in vitro tests: APS-HA samples (a, b, c, d); APS-HA/SPS-BG_Ca samples (d, e, f, g). The scale bar is the same for all images, secondary electrons.

Fig. 9 .

 9 Fig. 9. Evolution of the XRD spectra of the APS-HA samples (bold trace) and of those of the APS-HA/SPS-BG_Ca samples (thin trace) after in vitro tests. The spectra of the as-sprayed samples are included as a term of comparison. Symbols: *, larnite (Ca 2 SiO 4 ); ▲, TTCP; ▼, TCP; ■, CaO; □, CaCO 3 ; •, HCA. The peaks without symbol are due to HA.

Fig. 10 .

 10 Fig. 10. Evolution of the Raman spectra of the APS-HA samples (bold trace) and of those of the APS-HA/SPS-BG_Ca samples (thin trace) after in vitro tests. The spectra of the assprayed samples are included as a term of comparison.

Table 1

 1 Deposition parameters used to apply the pure HA coating by plasma spraying (APS-HA) and the bioactive glass topcoat by suspension plasma spraying (SPS-BG_Ca).

		APS-HA	SPS-BG_Ca
	Torch power [kW]:	21	40
	Gases
	Ar [slpm]:	48	45
	H 2 [slpm]:	2	7.5
	Spray distance [mm]:	100	70
	Torch speed [mm/s]:	500	750
	Scan step [mm]:	10	10
	Feedstock flow rate [g/min]:	10	30
	Injector diameter [mm]:	1	0.3
	Carrier gas flow rate [slpm]:	Ar, 3	-
	Static pressure [MPa]:	-	0.17
	Number of torch passes:	36	2 × 4 runs

A. Cattini et al. / Surface & Coatings Technology 240 (2014) 110-117

Acknowledgements

The Vinci Program (Italian-French University) support for Dr. A. Cattini is gratefully acknowledged.


The analysis of possible application of different working gases in cold gas spray method: dynamics and temperatures of particles

The paper analyzes the possibility of application of different working gases and their mixtures in cold gas spraying method (CGSM). The gases such as: (i) H2O vapors; (ii) CO2 ; and (iii) Ar were considered. On the other hand, the analysis included the gas mixtures of Ar, CO2 and N2 with H2O vapor. The temperature of analyzed gases was supposed to be T=800 K and its pressure P=2 MPa which are the typical values in CGSM. The calculation of the profiles along the convergent-divergent nozzle of gas and particle velocity was carried out with the use of Fluent software. The velocity profiles of all analyzed gases were compared to that of the standard gas in CGSM, nitrogen flowing in the same conditions. Finally, the velocities and temperatures of copper particles having diameters 20 and 50 μm in the tested gases flows were calculated and compared to that reached in nitrogen in similar conditions. The difference in temperatures and velocities were explained by sonic velocities, viscosities and thermal conductivities of selected gases.

Introduction

The cold gas spray method (CGSM) is one of the newest methods among thermal spraying techniques in which small unmelted powder particles are accelerated in the pressurised gas stream towards the substrate. The bonding is obtained due to excessive plastic deformation which occurs during particles impact. The bonding strength depends strongly on particles velocity which must be above so called critical velocity to enable deposition [START_REF] Schmidt | Development of a Generalized Parameter Window for Cold Spray Deposition[END_REF]. The rise of the particles velocity is possible by: (i) an increase of the temperature and pressure of the gas; (ii) a modification the nozzle geometry; and, (iii) proper choice of the gas type [START_REF] Champagne | The cold spray materials deposition process -Fundamentals and applications[END_REF]. The gases used in cold spray process should have high sonic velocity as well as high viscosity in order to transfer efficiently the momentum to powder particles. They should also be cheap, safe and easily available [START_REF] Champagne | The cold spray materials deposition process -Fundamentals and applications[END_REF]. Currently, nitrogen is very frequently used but when such materials as titanium are processed, the particles velocities do not reach the critical values. Consequently, more expensive gas such as helium must be used. The present paper analyses the possibility of application of gases able to replace helium namely steam, argon, carbon dioxide, as well as their mixtures. Water steam seems to be the most promising candidate as it has high sonic velocity and is cheap as easy accessible. The idea of using water vapor for accelerating particles was already mentioned by Nesterovich et al. [START_REF] Nesterovich | Process of Powder Compaction[END_REF] but needs more studies verifying usefulness of this gas.

Calculation model

The simulation of the gas flow and its interaction with particles was done using commercial software Fluent (ANSYS, Fluent Inc.,USA). The coupled implicit method was applied to simulate the gas flow. After having obtained the stable solution, the particles were introduced to the gas and their temperatures and velocities were calculating with the discrete phase model (DPM). The gas flow, its turbulent kinetic energy, and turbulent dissipation rate were modelled using second order upwind accuracy [START_REF] Ning | Numerical Study of In-flight Particle Parameters in Low-Pressure Cold Spray Process[END_REF].

The computational domain

The calculation axisymmetric area in 2D is shown in the Figure 1. The simulation was performed using convergent-divergent type of nozzle having cylindrical geometry with length of 255 mm and internal diameter at exit I.D.=4 mm. Argon, carbon dioxide, water steam and their mixture with pressure 2 MPa and temperature of 800 K were used as working gases.

The powder was fed axially. The atmosphere was simulated as cylinder of 30 mm diameter and 150 mm length starting from the nozzle exit.

Fig. 1 The computational domain

The boundary conditions

The stagnation pressure (absolute) at inlet was set to be P = 2.1 MPa and the temperature T = 800 K. For the exit area the pressure and temperature was 0.1 MPa and 293 K respectively. The nozzle was defined as smooth, not moving boundary without the heat flux with the walls. The gases were assumed to be compressible, viscous, turbulent and obeying the ideal gas laws. The latter is doubtful in the case of steam but as was proved that such assumption is acceptable and gives reasonable results [START_REF] Aphornratana | Theoretical and experimental investigation of a combine ejectorabsorption refrigerator[END_REF][START_REF] Sriveerakul | Performance prediction of steam ejector using computational fluid dynamics (I) Validation of the CFD results[END_REF]. The values of thermal conductivity, specific heat capacity were taken from the Fluent database and were assumed to be constant during simulation. The changes in dynamic viscosity with temperature were taken in account using Sutherland formula [START_REF]Ansys Fluent 14.5 User's Guide[END_REF]:

ߤൌߤ ቆ ܶ ܶ ቇ ଷȀଶ ܶ ܵ ܶܵ (1)
Where, Tref is reference temperature, ȝref -dynamic viscosity at the Tref, S -Sutherland temperature.

The turbulence of gas was modelled using standard k-H model which is computationally efficient and gives a fairly accurate prediction of the turbulence in the flow of powder particles in the gas stream flowing through the convergent-divergent nozzle [START_REF] Liao | Optimal design of a cold spray nozzle by numerical analysis of particle velocity and experimental validation with 316L stainless steel powder[END_REF][START_REF] Li | Optimal design of a novel cold spray gun nozzle at a limited space[END_REF][START_REF] Pattison | Standoff distance and bow shock phenomena in the cold spray process[END_REF].

The properties of gases mixtures were modelled using gas mixing laws. The molar ratio of gases in mixtures was 1:1.

The spherical copper particles with diameter 20 and 50 ȝm were used in modelling. The initial velocity was set to 5 m/s and temperature to 300 K.

Heat and momentum transfer between gas and particles

Particles were modelled as a discrete phase dispersed in a continuous phase and the motion equation was specified using the Lagrange formula [START_REF]Ansys Fluent 14.5 User's Guide[END_REF]. Excluding the gravity force, the balanced of forces acting on a particle is as follows:

ݑ݀  ݐ݀ ݉  ൌܨ ܦ ( 2
)
where FD -the drag force, mp -mass of particle, upparticle velocity.

The drag force can be expressed by the following relationship:

ܨ ܦ ൌ ͳͺߤ ߩ  ݀  ʹ ܥ ܦ ܴ݁  ʹͶ ൫ݑ ݃ െݑ  ൯ (3)
where ȝg -dynamic viscosity of the gas, ȡp -density of particle, dp -diameter of particle, CD -drag coefficient, ug -gas velocity, up -particle velocity, Rep is the Reynolds number:

ܴ݁  ൌ ߩ  ݀  หݑ ݃ െݑ  ห ߤ ݃ (4)
In case of spherical particle the drag coefficient was taken from Morsi and Alexander [START_REF] Morsi | An Investigation of Particle Trajectories in Two-Phase Flow Systems[END_REF]:

ܥ ൌܽ ଵ ܽ ଶ ܴ݁ ܽ ଷ ܴ݁ ଶ (5)
Where a1, a2, a3 are the coefficients applied over several ranges of Rep.

Since it was assumed that the temperature of particle is below the vaporization temperature and the temperature gradient inside the particle was negligible, the particle temperature during process may be taken as follows [START_REF] Ning | Numerical Study of In-flight Particle Parameters in Low-Pressure Cold Spray Process[END_REF][START_REF]Ansys Fluent 14.5 User's Guide[END_REF]:

݀ܶ ݐ݀ ൌ ݄ήܣ ήሺܶ െܶ ሻ ݉ ήܥ (6)
where: h -heat transfer coefficient of gas, Ap -outside surface area of a particle Te -gas temperature at the nozzle exit,Tp0 -temperature of powder before introduction to the nozzle, Cp -specific heat of particle.

Results and discussion

The velocity profile of water steam and that of sprayed particles are shown in Fig. 2. The analyzed copper particles reached on the outlet of the nozzle the velocities which ranged from 402 m/s to 571 m/s in steam and from 374 m/s to 512 m/s in nitrogen depending on their diameter. The both gases accelerated the particles to the greatest velocities and the calculation for argon and carbon dioxide in the same conditions gave the values of 350 m/s to 460 m/s and 341 m/s to 454 m/s respectively (Fig. 6).

Fig. 6 Comparison of particles velocities at the nozzle outlet for different working gases

The influence of working gas on particle velocity was much clearer for small particles. In case of particle with diameter 50 ȝm, the velocity in water steam was about 7% greater than in nitrogen. This difference increased to 12 % for 20 ȝm particles. After leaving the nozzle the particles velocities in water steam and nitrogen are retained and even increased slightly.

Outside the nozzle the shockwaves (oscillation of velocities) are visible for all gases. The sharp decline of the gas velocity follows the shockwave. The acceleration of particles takes place along the nozzle and outside of it. Thus, the highest particle velocity is obtained at the distance of 45 mm from the nozzle exit. The velocity rise is similar for all gasses and equal to about 20 m/s.

The temperature of all gases is similar in the convergent part of nozzle and decreases in its diverging part where the gases decompress. The decompression is the most intensive for argon. This results in the lowest temperature of gas, Tg = 428 K and of particles, Tp = 564 K to 627 K at the nozzle exit. The decompression is much less intensive for carbon dioxide. The temperature of this gas in the same position is Tg = 613 K and that of particles is Tp = 641 K to 709 K (Fig. 7). These temperatures are well below the melting point of copper.

Fig. 7 Comparison of particles temperatures for different working gases at the nozzle exit

The obtained results can be explained by the differences in following gases properties [START_REF] Dykhuizen | Gas Dynamic Principles of Cold Spray[END_REF]: (i) the ratio of specific heats, Ȗ (ii) the specific gas constant, R (the universal gas constant is divided by the gas molecular weight); and, (iii) dynamic viscosity, ȝ. The influence of these parameters can be explained using isentropic flow theory. Following this theory, the velocity of gas in the nozzle can be calculated from the following equation [START_REF] Dykhuizen | Gas Dynamic Principles of Cold Spray[END_REF]: ܸൌܯඥߛܴܶ

where M -Mach number, ߛ -specific heat ratio. As the Mach number depends only form nozzle geometry, by constant temperature, the increase of the velocity is possible only by using gas with higher specific heat ratio or specific gas constant.

Additionally the dynamic viscosity determines the drag force (see eq. 3-5).

Water steam has the highest specific gas constant and the particles sprayed with have, consequently, the highest gas velocity. On the other hand, smaller dynamic viscosity of this gas in comparison with nitrogen reduces considerably the momentum transfer to sprayed particles. The particles velocities in these two gases results, therefore, rather from the powder material density and particles size which influence the drag force acting on them as shown in [START_REF] Malachowska | The possible application of steam as working gas in cold gas spray method: analysis of momentum and heat transfers from gas to particles[END_REF]. Consequently, water steam is an efficient working gas to spray light materials having small size.

It should be mentioned that the condensation was not taken into account in used models. On the other hand, the steam would surely transform into water after leaving the nozzle and also on the nozzle walls (to a small degree). It is supposed that this effect would result in an increase of dynamic viscosity and, consequently, in higher particles velocities. However, the further studies are needed to take into account this phenomenon. The velocity of the powder is important because it influences the coatings' adhesion. Only the particles exceeding critical velocity, which is specific for each material, can be deposited [START_REF] Fukanuma | On-flight particle velocity measurements with DPV-2000 in cold spray[END_REF]. The critical velocity for copper on copper substrate was found to be 422 +/-45 m/s [START_REF] Raletz | Critical particle velocity under cold spray conditions[END_REF]. Consequently, all tested working gases proved to be sufficient to obtain copper coating.

The other important factor influencing the adhesion is particles temperature influencing its thermomechanical properties [START_REF] Ghelichi | Numerical simulation of cold spray coating[END_REF].

Many materials, such as metals and alloys, deforms easier in a plastic way with increasing temperature. The cold sprayed powders of such materials deforms at impact with substrate better and their bonding to substrate improves when particles temperatures increases. However this temperature may introduce oxidation what should be avoided when cold spray.

The temperature of the particles may be controlled by the shape of the nozzle by the selection of working gas. Thus, the use of argon, which has high ratio of specific heats, results in particles having much lower temperature than that at spraying with carbon dioxide. Thes calculated values for 20 ȝm copper particles are Tp = 627 K and 709 K respectively.

Fig. 7 Comparison of particles velocities for different gases mixtures

The effects th gas mixing on particles velocity is shown in Fig. 7. The calculated particles velocities were lower comparing to that obtained for pure steam. Therefore the pure water steam remains the most interesting potential candidate for working in CGSM.

The important residual problem which remains to be solved is the condensation of water on the substrate. One of the solutions is the strong overheating of steam and heating the substrate well above the steam condensation temperature.

Conclusions

The paper analyzes possibility of application of different gas for cold gas spray method in comparison to nitrogen, taken as a standard. The temperature and pressure were set as typical value used in this method. The commercial software Fluent was used to calculate the velocities and temperatures of particles at processing with different gases. The comparison shows that the highest velocities were obtained for pure steam. The highest temperature were reached by particles sprayed with carbon dioxide and the lowest -with argon. The further work will be carried out to estimate the different effects of steam condensation on the velocity, temperature of sprayed particles and the quality of obtained coatings.

  Fig 3 shows the similar profiles for nitrogen. The temperatures in Fig 4 and Fig 5 respectively.

Fig. 2 Fig. 3 Fig. 4 Fig. 5

 2345 Fig. 2 Profiles of gas and particles velocities in the nozzle axis for steam at 2 MPa and 800 K

Acknowledgments

Marcin Winnicki has contributed in the research of literature cited in this paper.