E Stratta

M Belliard
email: michel.belliard@cea.fr

THERMAL-HYDRAULIC STUDY OF PASSIVE SAFETY SYSTEMS BASED ON THE HYDRAULIC DIODE PRINCIPLE FOR THE MANAGEMENT OF LARGE-BREAK LOSS OF COOLANT ACCIDENTS

Keywords: Passive safety systems, In-vessel flow limiter, Advanced accumulator, Large-break LOCA

In the context of improvement studies for pressurized water reactors (PWR) in the post-Fukushima accident situation, passive safety systems have been identified as an interesting strategy for design-basis and severe accidents management. CATHARE code thermal-hydraulic studies, performed to evaluate the impact of two specific passive elements of safety systems during large-break loss-of-coolant-accidents (LB LOCA), are here presented. The two systems, namely the in-vessel flow limiter and the advanced accumulator, are based on the hydraulic diode principle: a singular pressure drop is added to the flow path by an induced fluid vortex, thus decreasing the mass flow. The in-vessel flow limiter is made up of some fins deflecting the fluid only in case of reversed flow, thus limiting the liquid mass exiting the primary circuit through the break. The advanced accumulator modifies its injection flow by switching the orientation of the injection pipe. After the reactor refill (high-rate injection phase), the flow regulation prevents the waste of the accumulator water: the lower the flow, the more lasting the injection. Parametric studies on the pressure-drop coefficient have been performed to evaluate each system's impact. The flow limiter coefficients validity domain has been verified by CFD studies (GENEPI code). Analyses on the maximum cladding temperature, the primary-circuit mass inventory and the core reflooding-front level have shown a benefit mainly depending on the pressure-drop coefficient value. The two systems' combined effect shows a great gain on the three mentioned parameters, even suggesting the possibility of reducing the active safety injection, thus reducing possible failures and investment costs.

INTRODUCTION

The present work sets in the context of Generation II and III Nuclear Power Plants innovation in the post-Fukushima accident situation. Since the French built reactors are mainly light water reactors (all the exploited ones anyway), the focus is set on these ones, notably the PWRs. Subsequently to the active safety systems failure after the tsunami in Fukushima, a great interest in passive safety systems arose. Passive systems being characterized by no need of external power forces, external input signals or moving mechanical parts, they have been detected as an appealing strategy in case of environmental disasters. Passive safety systems have gained (and are continuously obtaining) more importance in the last years, due to multiple factors. First of all there is the interest in a reactor which could handle LOCA accidents (still considered one of the most severe design basis accidents that may occur) without external help, under the condition of demonstrating the operation of those systems.

In this paper the focus is set on safety systems that work on the hydraulic diode principle. The hydraulic diode is a device that can be actuated under different forms. Basically the idea is the same as an electrical diode, where this component allows the passage of the electric current in one direction, but nearly blocks it in the opposite one. For the hydraulic diode the mechanism is the same but it is actuated by an increased pressure drop on one of the two sides of the flow orientation. This device has been applied to two different systems: the in-vessel flow limiter and the advanced accumulator. The standard PWRs do not employ any of these passive systems, but there are some examples of new Generation III reactors that adopt the advanced accumulator, notably the ATMEA1 reactor and the APR1400 (e.g. the South Korean reactor Shin-Kori3 and 4 under construction) [START_REF]Shin-Kori reactor characteristics[END_REF]. At the present day no designed reactor includes the flow-limiter between the vessel and the loops. Nevertheless, patents and studies have already been performed at the CEA in the past years and the interest towards this system is still alive. On each of these systems a parametric study has been done, to evaluate the efficiency of the device under LB LOCA conditions. LB LOCAs have been chosen as the reference accidental transients because of the rapid vessel emptying: these are the case in which it is the most important not to waste water. As that the main purpose of one of the two devicesthe in-vessel flow limiteris to reduce the amount of lost water, LB LOCA transients are the most suitable to highlight the system's impact. In this study we will consider a prototypal 3-loop PWR reactor (1150MWe, N4 SGs 1 ) that will be presented in Section 2 and modelled using the French thermal-hydraulic safety code CATHARE [START_REF] Barré | The CATHARE code strategy and assessment[END_REF]. In the two following sections we will describe specific hydraulic diode based safety systems, including their CATHARE model and accidental transient results on LB LOCA. Finally we will draw some conclusions and perspectives.

CATHARE REACTOR MODEL

CATHARE is a reference thermal-hydraulic safety code originally devoted to the study of different types of transients that can occur in water-cooled reactors (both standard operation and accidental transients). It is based on modules that are computed through a six-equation two-fluid model (mass, momentum and energy equations for the two phases) with the possibility of adding supplementary equations that model non-condensable gases and radiochemical elements. The reactor considered in the study is a PWR that takes inspiration from the ATMEA1 reactor [START_REF]Status Report for Advanced Nuclear Reactor Design -Report 99[END_REF], a jointventure of AREVA and MITSUBISHI. It has been modelled with the CATHARE code (cf. Figure 1), using axial, volume, boundary condition and specific modules (e.g. for the guillotine rupture).

Figure 1. CATHARE simplified reactor scheme view

Our vessel model includes the downcomer channel, the lower plenum, the core, the bypass and the upper plenum. The water enters the vessel in an annular zone that has been modelled through a small inlet volume and an axial downcomer, it is collected in the lower plenum and after that it ascends the core to cool the fuel elements: these last are 4.2 m long axial elements that model the average core and the hottest channel (point-kinetics neutronics modules constitute the fuel pins). A bypass element defines the amount of water that does not actually cool the fuel. The heated water is then gathered in a volume element that represents the upper plenum. From there, the water is then redirected towards the three loops. In order to take into account the parasite losses of water at the inlet of the vessel, a bypass towards the upper plenum has been defined. Specifically, on each intact looplocated just after the primary pumpsthere are three safety injection systems: a medium-head pump (MHSI priming pressure at 85 bar), a low-head pump (LHSI priming pressure at 20 bar) and an accumulator (initial pressure: 47 bar). The secondary circuit consists of three identical steam generators that are coupled with the three loops of the primary system. Each steam generator is made up of some axial modules representing the two specific parts of the downcomer and lower riser and the common upper riser. The steam dome volume element is connected by the steam line is to the steam head (just before the turbine). Some controllable BC modules model the ARE and ASG injections as well as the steam line output.

IN-VESSEL FLOW LIMITER

Description

Figure 2 illustrates the position and the design of the in-vessel flow limiter, which is applied on the connections between the hot leg of the loops and the vessel. The flow limiter is applied to fulfil two goals. First of all, it aims at reducing the water flow that exits the vessel in the direction of the break. Secondly, it is conceived to avoid the direct transfer of the intact loops injection water to the broken loop (and so eventually towards the break). From the intact cold leg, the water should flow into the downcomer and pass through the core to cool it before rise up and at last exit the circuit. In order to attain the first goal, some fins (composing the diode) are placed on the vessel-loops connections so that a large vortex is created in case of reversed flow. The pressure drop coefficient is so increased compared to the normal flow case (regular hydraulic path during normal operations). For the second objective there are three comparted slots that divide in three parts the annular downcomer (just one inlet/outlet per sector).

The CATHARE Model

To get a first insight of the hydraulic behaviour of the in-vessel flow limiter, our model consists in modelling the effect of the hydraulic diode through the variation of the pressure drop coefficient between the downcomer and the cold legs. In the standard CATHARE case, the geometric model is made up of three separate loops which converge onto a single vessel, connected by junctions onto which are defined the pressure drop coefficients. These coefficients are symmetric. This means that the code will compute the same value of pressure losses, the flow being direct or reverse. For the standard case the IN/OUT pressure loss coefficient is arbitrarily taken as k IN = k OUT = 0.78. In the first and simplest modification that has been applied to the data setwith the same geometry of the standard casethe pressure drop coefficients for the reverse flow (in opposition to the default flow direction, which is defined by the steady state flow direction) have been increased to simulate the diode.

CATHARE LB-LOCA transient analysis

A parametric study on the value of the pressure drop coefficient k OUT has been realised. CATHARE does not allow running CFD (Computational Fluid Dynamics) simulations including the effective geometry of the fins and the turbulence that derives from it, so it is really tough to evaluate the exact value to be set as coefficient. The input data for the parametric study will be assessed by CFD evaluations in order to demonstrate the feasibility of the device (cf. Section 3.4). The k OUT (coefficient for the reverse flow) has been incremented by applying various multiplicative factors: K = 5 -10 -15 -100, leading to the following pressure drop coefficients k OUT : 0.7783 -3.8915 -7.783 -11.6745 -77.83. Taking into account literature [START_REF] Gautier | Development of the safety systems for a simplified, low pressure, medium size PWR[END_REF], K =15 has been chosen as the reference, since it is the highest value that remains realistic (K = 100 is not likely to be feasible). Computations have been performed under two safety injection logics. The first case is run with only one accumulator per intact loop ("without MHSI/LHSI" case). The second case accounts for one MHSI pump and one LHSI pump per loop in addition to the accumulator ("with MHSI/LHSI" case). The three main effects related to the diode application have been identified into short-term and long-term effects. They will now be analysed in details, through some comparative graphs (cf. Figure 3, Figure 4 and Figure 5). The short term is referred to the first 30-50s of the transient. During this extent of time the diode helps to reduce the emptying of the primary circuit and the early-peak of the cladding temperature. Whichever the safety injection case, the results of the first 30s are the same in the two transients (with or without SI pumps), because the additional injection sources begin to work with a small delay (at 32s), thus not influencing the first part of the accident. In Figure 3-left it is possible to see the first major effect of the diode. Comparing only the case with K = 15 to the standard case (K = 1), the minimum value of the primary mass is clearly higher in the case where the diode is in place. During the first 30s the curves representing the behaviours "without MHSI/LHSI" and "with MHSI/LHSI" are superposed and just at the pump activation moment there is a consistent difference between the two cases. Figure 3-right shows the core inlet flowrate: the second beneficial effect of the hydraulic diode is notably the possibility to have a non-negligible flowrate that enters the core during the first instants of the transient. The red curve (diode with K=15) well describes this behaviour if compared to the standard case. Figure 4 shows the temperature profiles of the hot pin cladding ("with MHSI/LHSI" case), with a zoom on the first 100s, detailing the short term physical impact of the diode on the system. The cladding temperature peak is consistently decreasing with respect to the increase of the multiplicative coefficient: the blue curve (K=15) experiences a peak reduction of almost 100°C compared to the standard case. This means that, if a diode with a k OUT = 11,6745 can be created, the gain in terms of temperature decrease will be appreciable. After the peak, a consistent decrease of the temperature is due to a reverse flow that comes from the top of the core. The cladding cools down and the system begins to store the injected water in the lower plenum. It is for this reason that after ~10s the temperature rises again, even with the accumulators' injection. Until the lower plenum and the downcomer are not full, it is not possible to cool the core. This is achieved at about 50s, and finally the water can cool the fuel elements. On the long term the transient differs depending on the safety injection logic: in this section we focus on the case "with MHSI/LHSI" (cf. Figure 4 and Figure 5), the detailed results concerning the case "without MHSI/LHSI" can be found in Section 4.3. The core re-flooding is the most important parameter to be analysed, from which also the cladding temperature is dependent. If the re-flooding is attained, the meltdown of the core will be prevented (even if maybe just temporarily), as well as all the consequences that follow. Figure 5 shows that the hydraulic diode improves the core reflooding: thanks to the flow limiter, the top of the core is attained earlier by the quench front. Compared to the standard case, the blue curve (K=15) anticipates nearly of 40s, which is absolutely not negligible during an accident in a nuclear reactor. Moreover, when the water level is high enough to completely cool the cladding, the cladding temperature drops (cf. Figure 4). In Section 4.3, the figures of the "without MHSILHSI" case show that clearly the liquid water does not reach the total height of the fuel. However, thanks to the diode, the attained level is higher than the standard case. As a consequence, the limited water safety injection cannot avoid (on the long term) the temperature divergence under the heating of the decaying residual heat, even if the effect of the diode is present also on the long term: the heating is a bit delayed in time and lower than in the standard case. The obtained results are coherent with the literature ones [START_REF] Gautier | Development of the safety systems for a simplified, low pressure, medium size PWR[END_REF], even if the comparison is not always easy because the operation pressure of the reactors is different, the safety injection is not always the same and finally because the absolute values for the pressure drop coefficients are not known. We can notice that, due to nonlinear-behaviour of the final reflooding level and consequently of the final temperature, the choice of the multiplicative coefficient K is to be optimized. Anyway, it is clear that the diode effect is always beneficial compared to the standard case. Table I presents results for a diode that increases the pressure drop coefficient 15 times. In case of full SI, the re-flooding is always reached because of the great injection of water due to the pumps. Anyway, the increase of the pressure drop coefficient k OUT allows the maximisation of the primary mass (almost 3 tons more) and the reduction of the reflooding delay (about 40s). On the other hand, it is remarkable that, even with two accumulators only, the diode-like device allows to reduce the short term cladding temperature peak (almost 100°C) and avoid the core melt-down. In terms of safety it is clearly an interesting benefit.

Table I. In-vessel flow limiter main results

Standard

(GENEPI) CFD Pressure drop coefficients Assessment

The motivation of CFD studies concerning the in-vessel flow limiter is the consolidation of the CATHARE results based on parametric studies of the pressure-drop coefficient k OUT . Section 3.3 highlights the benefice brought by the value k OUT = 11.6745 (or multiplicative coefficient K = 15). As the value of the pressure-drop coefficient is linked to the fan design geometry, we have to consider an optimized design to assess the pressure-drop magnitude values used in the CATHARE study. Optimizing the in-vessel limiter design needs a huge number of computations needing small CPU times. For this reason, we choose to use the GENEPI2 code [START_REF] Grandotto | Steam generator gwo-phase-flow numerical simulation with liquid and gas momentum equations[END_REF]. It was designed for the steam-generator two-phase flow steady-state 3D computations through the resolution of the mass and momentum balance equations of a dilatable liquid/steam mixture. Moreover, it incorporates the possibility to model thin no-penetration obstacles using Immersed Boundary Conditions (IBCs [START_REF] Belliard | Fictitious domain methods for two-phase flow energy balance computations in nuclear components[END_REF]): we consider simulations over a full computational domain including the in-vessel flow limiter and re-introduce its presence adding local external forces (fictitious domain method) in order to get a fast estimation of the pressure-drop coefficient.

The goal of this section is to get some preliminary elements to answer to the question: is the value k OUT = 11.6745 realistic? In this preliminary study, we restrict ourselves to mono-phasic liquid steadystate flows and to an approximated (not optimized) in-vessel limiter design. We only model one third of the down-comer, including one cold-leg entry in the vessel and one hot-leg pipe. The CFD computational domain is based on the CATHARE space discretization, cf. Figure 6. It is a simplified planar geometry to which was added the broken cold-leg nozzle starting 0.93 m before the down-comer. It extends up to 1.8974 m below the cold-leg axis. The geometrical dimensions are taken from the CATHARE data set, except for the thickness of the planar part. A compromise between the volumes and the radius leads us to set this data to about 0.20 m. For the purpose of a mesh convergence study, three meshes M 1 , M 2 and M 3 were built involving N 1 = 6.080, N 2 = 48.640 and N 3 = 164.160 elements respectively. Each mesh i can be characterized by a space-step index (N i ) -1/3 . The ratio between two consecutive space-step indexes is 2.0 (M 1 →M 2 ) and 1.5 (M 2 →M 3 ).

𝐾 𝑔𝑙 = 2(〈𝑃 𝑖𝑛 〉 -〈𝑃 𝑜𝑢𝑡 〉)/(〈𝜌 𝑐𝑙 〉〈|𝑉 𝑐𝑙 𝑛 𝑐𝑙 |〉 2 ) (1)
where the symbol <•> denotes an area average, P in the pressure at the in-flow surface and P out the pressure at the out-flow one. We compare this global pressure-drop coefficient to the CATHARE estimation of the same quantity K gl,cath . From a discrete-space point of view, considering the three meshes M i = M 1 …M 3 , we denote by K Mi the computed approximations of K gl .

As the GENEPI2's turbulence viscosity model 𝜇 T is quite rough (the Schlichting model [START_REF] Grandotto | Steam generator gwo-phase-flow numerical simulation with liquid and gas momentum equations[END_REF]: 𝜇 T = a S |G| L T with G the mass flux), a parametric study is performed on the coefficient a S and the turbulence characteristic length L T . The characteristic length is related to the biggest eddy structures. As reference, we choose L T ~ 1 m (azimuth scale in the down-comer) for the nominal-operation flow direction and L T ~ 0.3 m (radial scale in the down-comer) for reverse flow direction toward the broken cold leg. In GENEPI2, the standard value for the Schlichting coefficient is a S = 0.047. We present some field distributions on the finest mesh concerning the mixture pressure and the mixture velocity for the two studied flow configurations: the nominal-condition one and the reverse-condition one (LB LOCA). . This must be confronted with the primary-pump characteristics to check if this overhead can be offset.  However, for reverse flows, the global pressure-drop coefficient is multiplied by almost a factor two in case of flow limiter: [3.6; 3.7] → [5.4; 5.6]. This is a bit less that the increase computed by CATHARE in the case k OUT = 11.6745: a factor three (1.3 → 3.9). But the CFD simulation reaches pressure-drop coefficients bigger than those needed by CATHARE to get a beneficial effect on the LB LOCA transient: about 5.5 to compare to 3.9. Then, we can conclude that the CATHARE flow-limiter pressure-drop coefficient k OUT = 11.6745 is consistent with those issued from our preliminary GENEPI2 study. Nevertheless, considering the limitations related to the turbulence description, this conclusion needs to be consolidated by body-fitted CFD studies with more precise turbulence models.

ADVANCED ACCUMULATOR

Description

The second passive system using a hydraulic diode is the "advanced" accumulator. Figure 9 and Figure 10 show a simplified scheme of the design and the mode of operation. The special feature of this accumulator is the double-regime injection flowrate. At first the flowrate is huge, allowing the refill of the circuit, but after some time, the flowrate is automatically switched to a lower value, just to compensate the water that exits the reactor through the break.

Figure 9. Advanced accumulator [10]

The actuation mode is passive, because it does not need any mechanical activation. When the water level is higher than the height of the vertical standpipe (cf. Figure 9), there is a big flowrate that descends into the flow damper and is then injected directly in the circuit. When, instead, the level is lower, the injection is produced only by the horizontal small flowrate pipe that is placed on the limit of the flow damper. In this way, this particular injection point will induce a vortex in the flow damper, thus defining the "vortex chamber". As in the flow-limiter situation, the vortex will increase the pressure drops and so reduce the injected flow.

Figure 10. Advanced accumulator injection law [10]

The interest of this advanced accumulator lies within the possibility of reducing the waste of water. Normally the safety injection of an accumulator is done as in the left graph of Figure 10, which means that the water is rapidly injected and in this way early consumed. For this reason the low-pressure injection pumps must start very soon in the evolution of the accident. The advanced accumulator injection follows instead the curve of the graph on the right of the same figure: changing the regime, the remaining amount of water is spent "smartly", thus allowing a delayed time of the pump employment.

CATHARE Model

The CATHARE model consists of a standard accumulator which changes the pressure drop coefficient on the injection line outlet. For the large flow regime model, the standard accumulator pressure drop coefficient value (3.813 i.e. as in the computations of Section 3.3) has been used. For the small flow regime model, this coefficient is changed to a lower one as soon as the prescribed low level is reached. In the transient block some flag indexes have been defined to test the water level inside the accumulator. This is done through the evaluation of the total mass M acc that remains and the hypothesis that the fluid density ρ water (and the accumulator base section S acc ) is constant:

𝐻 = 𝑀 𝑎𝑐𝑐 /(𝜌 𝑤𝑎𝑡𝑒𝑟 𝑆 𝑎𝑐𝑐 ) (2)
In order to appropriately define this system, the two degrees of freedom (H of the standpipe and the new pressure drop coefficient K sing ) have to be fixed. The regime switch level should be the one corresponding to the injected volume that allows the complete filling of the lower plenum and the whole downcomer. It is linked to the fact that after the complete filling of the lower plenum and the downcomer, the water that is sent into the circuit exits directly from the broken loop, thus being wasted.

The switch level H can be estimated by at least two approaches, one based on post-transient analysis (it can give the downcomer and lower plenum complete filling time, the remaining water volume in the accumulators and, hence, the level H) or on a priori geometrical considerations (for instance the sum of the geometric volumes of downcomer and lower plenum).

The pressure drop coefficient corresponds to an injected flowrate which maintains a natural circulation condition, by compensating the core vaporised steam flowrate. This allows us to design the pressure drop coefficient K sing by recalling the expression of the singular pressure drops:

|∆𝑃| = (𝐾 𝑠𝑖𝑛𝑔 𝐺 2 )/2𝜌 (3)
We consider the density ρ as constant (the temperature and pressure in the accumulator do not change in time). For the same pressure drop, it is necessary to increase the K sing so that the flowrate G decreases.

There is a quadratic inverse proportion that links the two variables G and K (where K 1 and G 1 are referred to the standard accumulator):

𝐾 1 𝐺 1 2 = 𝐾 𝑠𝑖𝑛𝑔 𝐺 2 2 → 𝐾 𝑠𝑖𝑛𝑔 = 𝐾 1 𝐺 1 2 /𝐺 2 2 (4)
The value K sing = 236 corresponds to a flowrate G 2 = 100 kg/s. It will be the reference coefficient for this study.

CATHARE LB-LOCA analysis

We consider here the same LB LOCA transient as in Section 3.3 and we analyse the case "without MHSI/LHSI". We retain the switch level based on previous computations post-transient analysis.

We are interested in evaluating the performances of the reactor equipped with the in-vessel flow limiter (K = 15) and the advanced accumulator (K sing = 236). In order to highlight the interest of the two allied systems, we provide comparison to the cases with standard or advanced accumulators (K sing = 236) without flow-limiter on the vessel (cf. Figure 11 on the left column).

After the presentation of the case with K sing = 236, as experimental data or CFD simulations are missing, we will present a parametric study concerning this pressure drop coefficient (cf. Figure 11 on the right).

On the following graphs are compared the ACCU flowrate, the quench front level and the cladding temperature. The best possible case contemplates the conjoint effect of the flow-limiter and the advanced accumulator (cf. Table II). First of all the accumulator injection time is more than doubled (130s = +135%), which means that the pumps can start with a bigger delay from the beginning of the accident.

Secondly, the reflooding level is increased of almost 1m (+34%), even if the top of the fuel is not completely covered. Finally, the gains on the cladding temperatures are perceptible: the early peak is reduced of about 100°C (-10%) and almost 600°C (-43%) on the long term.

Concerning the parametric study, the standard accumulator case (K sing = 3.813) has been compared to the following ones: K sing = 236, 120 (50% of K sing_ref ) and 60 (25% of K sing_ref ). As a whole, the results are slightly dependant on the pressure drop coefficient of the accumulator. Anyway the cladding temperature histories show an interesting point. While at the end of the transient the temperature is clearly lower for the case with the lowest injected flowrate, during the middle part of the accident the situation is inversed. This is due to the fact that the reduced flowrate is not strong enough to rapidly cool the core. It will take more time. Actually this is not a problem, because the temperature is anyway lower than the earlier experienced peak and this evolution leads to a final temperature that is more interesting for safety reasons. Table III sums up the results concerning the three most important features of this study: the accumulator emptying time, the core reflooding level and the hot pin cladding temperature at the end of the transient.

As it is possible to point out, the advanced accumulator is always beneficial in case of accident and the effect is proportional to whichever coefficient in the range [0, 236]. Even considering the smallest coefficient (K sing = 60, whose feasibility still needs verification), it allows a +5% accumulator emptying delay (and clearly a bigger delay before the need of safety injection pumps), it helps reflooding (+10%) and finally it reduces the cladding temperature (-20%), leaving a higher margin before the melt down.

CONCLUSIONS AND PERSPECTIVES

The study presented in this paper has been focused on the feasibility analysis of two passive safety systems that employ hydraulic diode devices: the in-vessel flow-limiter and the advanced accumulator. Considering a prototypal PWR, we have investigated the impact of these two passive devices on LB LOCA transients, using the CATHARE code. We have presented the CATHARE model of our prototypal PWR, of the in-vessel flow-limiter and the advanced accumulator. As a whole, each safety component has a beneficial impact on the transient, which is augmented when conjointly used. The benefices are to be found in term of reduction of peak cladding temperature, reflooding time and level. On the short term, the flow limiter helps the cladding temperature peak reduction, while the benefice on the long term depends both on the limiter (-22% if taken alone) and the accumulator (-30% if taken alone). However, the presented results have to be consolidated to qualify the physical achievability and economic feasibility of these passive safety components. Even if the present results are strongly encouraging, especially the conjoint use of the in-vessel flow limiter and the advanced accumulator, further CFD studies are needed to assess the input parameters of the study (in-vessel flow limiter and advanced accumulator pressure drop coefficients,…), by determining the best attainable geometrical designs. Let's notice in particular that the results on the flow-limiter are overall consistent with the mentioned literature [START_REF] Gautier | Development of the safety systems for a simplified, low pressure, medium size PWR[END_REF]. The two main effects that had already been enlightened by this literature are the presence of a core flowrate in the early instants of the transient and the consequent peak reduction for the cladding temperature. Actually, the studies are not easily comparable because of the different assumptions that are made at the beginning of the analysis: different primary pressures, dissimilar description of the safety injection systems and unknown absolute values of the pressure drop coefficients. Anyway, the resulting tendencies are in accord with the CEA previous examinations: for example, the peak reduction is weaker, but the overall evolution of the temperature is the same. Besides all the technical results that derive from the CATHARE computations, it is interesting to investigate the perspectives in term of active safety system reduction, which reduce investment costs and reactor operation, leaving unchanged (or even improved) the safety level. For example, analysing the obtained results in case of LB LOCA, even in the case "without MHSI/LHSI", the combined effect of the flow-limiter and the advanced accumulator, almost allows reflooding. Hence, the reduced requirement on the pumps efficacy (mass flowrate, …), forecasts a design that eliminates the large flowrate LP sources.

Figure 2 .

 2 Figure 2. In-vessel flow limiter position and design[START_REF] Belliard | Fictitious domain methods for two-phase flow energy balance computations in nuclear components[END_REF]

Figure 3

 3 Figure 3 On the left: Primary mass reduction -On the right: Core inlet flowrate

Figure 4 .Figure 5 .

 45 Figure 4. Parametric cladding temperature evolution with zoom -"with MHSI/LHSI" case

Figure 6 .Figure 7 .

 67 Figure 6. Part of the PWR geometry taken into account by the CFD study

Figure 7 .

 7 For this study, the pressure range considered in GENEPI2 is only bar. It differs from the pressure addressed in the CATHARE simulations of the nominal operation or of the LOCA, but it is the pressure drop that is important. Mass flux is imposed on the in-flow surface and pressure on the out-flow surface. The in-flow and out-flow surfaces are the down-comer and the cold-leg nozzle surfaces depending on the main flow direction (default, as in normal operation, or reverse, as in LB LOCA situation). These values are given by the CATHARE computations. For the nominal operation (default flow direction), the imposed mass flow rate equals 4.690 kg/s (cold leg) and the pressure equals 50 bar (down-comer). For the LB LOCA (reverse flow direction), the imposed mass flow rate equals 5.200 kg/s (down-comer) and the pressure equals 50 bar (cold leg). On the wall, we consider a slip wall boundary condition. Concerning the used methodology to compare the CATHARE results to the CFD ones, we do not compare the value of the local CATHARE pressure-drop coefficients k OUT to CFD estimation. But, we consider the global pressure-drop coefficient K gl computed between the in-flow/out-flow sections of the CFD computational domain marked in red in Figure7. Taking as reference the normal mass-flux at the beginning of the cold leg nozzle (G cl • n cl = ρ cl V cl • n cl with n cl the surface normal, V cl the flow velocity and ρ cl the flow density), we define K gl as:

Figure 8 .

 8 Figure 8. Pressure and velocity with flow limiter in normal operation (top) and LB LOCA (bottom)

Figure 11 .

 11 Figure 11. Accumulator injection flowrate -Reflooding level -Hot pin cladding temperature On the left: combined effect of the in-vessel flow limiter and the advanced accumulator On the right: parametric study for the accumulator pressure drop coefficient

Table II . Passive safety systems results

 II

		Standard ACCU	Advanced ACCU 236	Advanced ACCU 236 + Flow-limiter
	Time end ACCU [s]	55	126	130
	Reflooding level [m]	2.80	3.23	3.74
	Tclad @short-time	916	916	820
	Tclad @600s [°C]	1360	958	778

Table III . Advanced accumulator results (with in-vessel flow limiter K=15)

 III

		Standard	Advanced	Advanced	Advanced
		ACCU	ACCU 60	ACCU 120	ACCU 236
	Time end ACCU [s]	60	90	106	130
	Reflooding level [m]	3.22	3.57	3.67	3.74
	Tclad @600s [°C]	1063	869	833


PhD) a , R. Jennane T Janvier

MD, PhD Eric Lespessailles
email: eric.lespessailles@chr-orleans.fr

Subchondral Tibial Bone Texture Predicts the Incidence of Radiographic Knee Osteoarthritis: Data from the Osteoarthritis Initiative

Keywords: fractal analysis, trabecular bone texture, knee osteoarthritis, subchondral bone, radiographs

Objectives: To evaluate whether trabecular bone texture (TBT) parameters measured on computed radiographs could predict the onset of radiographic knee osteoarthritis (OA).

Materials and Methods:

Subjects from the Osteoarthritis Initiative with no sign of radiographic OA at baseline were included. Cases that developed either a global radiographic OA defined by the Kellgren-Lawrence (KL) scale, a joint space narrowing (JSN) or tibial osteophytes (TOS) were compared with the controls with no changes after 48 months of follow-up. Baseline bilateral fixed flexion computed radiographs were analyzed using a

Introduction

Although osteoarthritis (OA) is the most common disorder of the musculoskeletal system, there is still no cure and treatment options are limited, leading to an increasing number of total joint replacements [START_REF] Bitton | The economic burden of osteoarthritis[END_REF][START_REF] Losina | Lifetime medical costs of knee osteoarthritis management in the United States: Impact of extending indications for total knee arthroplasty[END_REF] worldwide. The underlying mechanisms through which this debilitating disease 1 occurs and progresses have not been fully elucidated yet [START_REF] Loeser | Osteoarthritis: A disease of the joint as an organ[END_REF]4 . Although the progression of the disease has been widely studied using several modalities [START_REF] Kraus | Subchondral bone trabecular integrity predicts and changes concurrently with radiographic and magnetic resonance imaging-determined knee osteoarthritis progression[END_REF][START_REF] Bruyere | Subchondral tibial bone mineral density predicts future joint space narrowing at the medial femoro-tibial compartment in patients with knee osteoarthritis[END_REF][START_REF] Roemer | Osteoarthritis Year in Review 2014: Imaging[END_REF][START_REF] Mackay | Subchondral bone in osteoarthritis: Association between MRI texture analysis and histomorphometry[END_REF][START_REF] Janvier | Subchondral tibial bone texture analysis predicts knee osteoarthritis progression: data from the Osteoarthritis Initiative[END_REF] , the genesis of the pathology has been much less investigated [START_REF] Liebl | Early T2 changes predict onset of radiographic knee osteoarthritis: data from the osteoarthritis initiative[END_REF][START_REF] Shamir | Early detection of radiographic knee osteoarthritis using computer-aided analysis[END_REF] .

Considering the growing body of evidence that showed the role of the subchondral bone [START_REF] Loeser | Osteoarthritis: A disease of the joint as an organ[END_REF][START_REF] Roemer | Osteoarthritis Year in Review 2014: Imaging[END_REF][START_REF] Kwan Tat | Targeting subchondral bone for treating osteoarthritis: what is the evidence? Best[END_REF][START_REF] Burr | Bone remodelling in osteoarthritis[END_REF] and the concept of crosstalk between articular cartilage and subchondral bone tissue [START_REF] Burr | Bone remodelling in osteoarthritis[END_REF] , we investigated the ability to predict the onset of knee OA based on observation of the subchondral bone. Trabecular bone texture (TBT) has been used for a long time to characterize altered trabecular bone [START_REF] Caldwell | Fractal dimension as a measure of altered trabecular bone in experimental inflammatory arthritis[END_REF] and is recognized to be a promising imaging biomarker [START_REF] Kraus | Trabecular morphometry by fractal signature analysis is a novel marker of osteoarthritis progression[END_REF] . Moreover, TBT parameters (TBTp) including fractal descriptors have been consistently reported to be able to predict the progression of radiographic knee OA [START_REF] Kraus | Subchondral bone trabecular integrity predicts and changes concurrently with radiographic and magnetic resonance imaging-determined knee osteoarthritis progression[END_REF][START_REF] Kraus | Trabecular morphometry by fractal signature analysis is a novel marker of osteoarthritis progression[END_REF][16][17] but also to be associated with the risk of total knee joint replacement 18 . Although conventional radiography has been considered unable to identify early relevant cartilage and bone changes associated with the onset of OA, TBT analysis could provide a second life to digital radiography and should provide new interesting information on the pathogenesis and occurrence of knee OA development 19 .

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Therefore, in the present longitudinal study conducted over 48 months of participants in the Osteoarthritis Initiative (OAI) we analyzed the predictive ability of TBTp measured on computed radiographs to predict the incidence of radiographic knee OA. There are very few studies [START_REF] Shamir | Early detection of radiographic knee osteoarthritis using computer-aided analysis[END_REF]19,20 that have investigated the genesis of radiographic knee OA through the TBT, yet, they showed that texture parameters combined with machine learning could predict knee OA initiation. Our contribution resides in the confirmation of these results on a large and well phenotyped cohort (OAI), using a reduced number of features reflecting the bone topology and considering a short-term outcome. A statistical approach robust to overfitting using cross-validation and a parsimonious reduction of the models was applied. Contrary to 19 , subjects with doubtful OA (KL=1) were excluded in the present work in order to focus on healthy at risk patients. In [START_REF] Shamir | Early detection of radiographic knee osteoarthritis using computer-aided analysis[END_REF]20 , they considered a long range outcome (20 years follow-up) while this study focuses on a short range outcome with a 4-year prediction. This work intends to identify patients with an elevated short-term risk to develop radiographic knee OA.

Material & Methods

Patients

The data used for the present study were obtained from the OAI database [START_REF] Nevitt | The Osteoarthritis Initiative: Protocol for the Cohort Study[END_REF] , which is available for public access at http://www.oai.ucsf.edu/. This database contains bilateral X-ray fixed flexion knee radiographs from 4796 patients followed since 2008. Patients were

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
5 recruited in 5 different centers and had annual bilateral fixed flexion knee radiographs using different equipment and modalities among the centers but all following the same acquisition protocol using the fixed flexion positioner SynaFlexer [START_REF] Kothari | Fixed-flexion radiography of the knee provides reproducible joint space width measurements in osteoarthritis[END_REF] . For this study, the radiographic scoring used was provided by the Boston University X-ray reading center [START_REF]Central Reading of Knee X-Rays for Kellgren & Lawrence Grade and Individual Radiographic Features of Tibiofemoral Knee OA[END_REF] . This assessment dataset includes both the Kellgren-Lawrence (KL) [START_REF] Kellgren | Radiological Assessment of[END_REF] and the Osteoarthritis Research Society International (OARSI) [START_REF] Altman | Atlas of individual radiographic features in osteoarthritis, revised[END_REF] semi-quantitative scores for radiographs from baseline to month 48, limiting the analyses to the 3421 graded knees in 1789 patients. As the present work focused on the initiation of knee OA, we excluded the prevalent doubtful OA knees (e.g KL = 1 at baseline), resulting in a subset of 514 knees from 471 patients. No standard protocols were provided for digitization of the film radiographs, resulting in considerable variability among the digitized images (variations in resolution, dynamic range, field of view, alignment, etc.) and they were therefore excluded.

The final subset included 344 knees in 319 patients. Only computed radiographs (CR) were analyzed. Digital radiographs (DR) were not considered in this study as none was available in the OAI at baseline. The complete pruning process is detailed in Figure 1 and the resulting population characteristics are detailed in Table 1.

OA grouping

Herein, the diagnosis of radiographic knee OA was defined as a KL score ≥ 1 according to the OAI definition [START_REF] Nevitt | The Osteoarthritis Initiative: Protocol for the Cohort Study[END_REF] . The onset of OA was investigated in patients with initial KL grade 0

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

(i.e. both joint space narrowing (JSN) and tibial osteophyte (TOS) grade 0) using three definitions based on these scores and described in Table 1. Knees with a KL going from 0 at baseline to 1 at 48 months (16 samples) were included as controls in order not to bias the results compared to a real case scenario. However, the same computations were achieved after removing these 16 samples, confirming the results presented here with non-significantly higher predictive values (data not shown).

Processing of radiographs

Radiographs were provided by the OAI as anonymous 16-bit DICOM images from which information such as the modality and the resolution were extracted. Each radiograph (right and left knees) was evaluated independently. Based on the results obtained in a previous study on the impact of ROI on OA characterization [START_REF] Janvier | ROI impact on the characterization of knee osteoarthritis using fractal analysis[END_REF] we chose to use the whole tibial bone area. The ROI extraction was performed by the principal investigator of the study, T. Janvier, as detailed in our previous study [START_REF] Janvier | Subchondral tibial bone texture analysis predicts knee osteoarthritis progression: data from the Osteoarthritis Initiative[END_REF] . The resulting patchwork of 16 ROIs is presented in Figure 2.

Each ROI side length was proportional to the tibia width to avoid misplacement due to anatomical variations in the knee widths. In our subset, radiographs presented different pixel spacing ሼ150; 158; 169; 170; 200ሽ µm without significant differences between groups, giving ROI side lengths ranging from 9 to 11 mm.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 7

Trabecular bone texture analysis

Fractal analysis has been regularly used in the past for TBT characterization [START_REF] Kraus | Trabecular morphometry by fractal signature analysis is a novel marker of osteoarthritis progression[END_REF]18,[START_REF] Messent | Fractal analysis of trabecular bone in knee osteoarthritis (OA) is a more sensitive marker of disease status than bone mineral density (BMD)[END_REF][START_REF] Wolski | Directional fractal signature methods for trabecular bone texture in hand radiographs: Data from the Osteoarthritis Initiative[END_REF][START_REF] Buckland-Wright | Subchondral bone changes in hand and knee osteoarthritis detected by radiography[END_REF] . It has been shown to be linked to the 3D microarchitecture of the bone [START_REF] Pothuaud | Fractal Dimension of Trabecular Bone Projection Texture Is Related to Three-Dimensional Microarchitecture[END_REF][START_REF] Jennane | Estimation of the 3D self-similarity parameter of trabecular bone from its 2D projection[END_REF][START_REF] Jennane | Fractal analysis of bone X-ray tomographic microscopy projections[END_REF][START_REF] Hirvasniemi | Correlation of Subchondral Bone Density and Structure from Plain Radiographs with Micro Computed Tomography Ex Vivo[END_REF] and to be robust to different common artifacts of X-ray imaging such as radiographic magnification, projection [START_REF] Lynch | A robust and accurate method for calculating the fractal signature of texture in macroradiographs of osteoarthritic knees[END_REF][START_REF] Winzenrieth | Three-Dimensional (3D) microarchitecture correlations with 2d projection image gray-level variations assessed by trabecular bone score using high-resolution computed tomographic acquisitions: Effects of resolution and noise[END_REF] and inhomogeneous contrast [START_REF] Pothuaud | Fractal Dimension of Trabecular Bone Projection Texture Is Related to Three-Dimensional Microarchitecture[END_REF][START_REF] Lynch | A robust and accurate method for calculating the fractal signature of texture in macroradiographs of osteoarthritic knees[END_REF] . For these reasons, images were processed directly from the DICOM files without any preprocessing.

In the present study, the fractal parameter (H) was computed using the variogram [START_REF] Bousson | Trabecular bone score (TBS): available knowledge, clinical relevance, and future prospects[END_REF] , which gives a measure of how much two pixels from a radiograph will vary in intensity depending on the distance between them. The H parameter is computed as the slope of the variogram in a log-log scale and quantifies the intuitive notion of the roughness of a texture (i.e. complexity of the texture evolution with the scale of observation).

The piecewise fractal nature of the TBT reported in [START_REF] Harrar | Piecewise Whittle estimator for trabecular bone radiograph characterization[END_REF][START_REF] Perrin | Piecewise fractional Brownian motion[END_REF] and confirmed by the variogram shapes in Figure 3 was also taken into account. The cut-off scale was observed around 500 µm on the empirical variograms and two fractal parameters were extracted: ܪ ஜ and ܪ corresponding to the texture complexity computed for the two micro (µ-scale) and the milli (m-scale) scales of observation respectively under 400 µm and above 600 µm. These two parameters ሼܪ ஜ ; ܪ ሽ were computed along both the vertical (90°) and horizontal (0°) axes corresponding to the main loading compression and tension directions, respectively. Finally, TBTp represent a vector of 4 descriptors ሼܪ ஜ|°; ܪ ஜ|ଽ°; ܪ |°; ܪ |ଽ°ሽ computed on 16 ROIs (Figure 2), resulting in 64 scalar parameters.

Statistical analysis

We compared the clinical covariates (age, gender, BMI) in patients with and without incident radiographic knee OA using the Mann-Whitney and Chi-square tests for numerical and categorical factors, respectively.

The predictive ability of TBTp was evaluated using a logistic regression through different models including combinations of the clinical covariates and the TBTp. To avoid the bias induced by a large number of parameters, these models were reduced using an automatic iterative stepwise optimization algorithm based on the Akaike Information Criterion (age, gender and BMI were not forced to be included in the final models). To avoid overfitting effects, the models used were evaluated using a 10-fold cross-validation repeated 300 times.

The performance of each model was evaluated through the receiver-operating characteristic (ROC) curves using the area under the curve (AUC) as a global performance criterion. The diagnostic odd ratio (DOR) was also computed [START_REF] Glas | The diagnostic odds ratio: a single indicator of test performance[END_REF] to investigate the relevance of the different models. DOR provides the ratio between the odds of a case being well classified against the odds of a control being wrongly classified as a case. It focuses on patients predicted as initiators and indicates if they were indeed true positives (DOR > 1) or false positives (DOR

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 9 < 1)
. This indicator is used to characterize the relative risk of actual radiographic knee OA incidence in a subpopulation predicted as initiators.

All statistical analyses were performed using R 41 (version 3.1.2) and the packages MASS [START_REF] Venables | Modern Applied Statistics with S. Fourth[END_REF] (version 7.3), Caret [START_REF] Kuhn | caret: Classification and Regression Training[END_REF] (version 6.0) and pROC [START_REF] Robin | pROC: an open-source package for R and S+ to analyze and compare ROC curves[END_REF] (version 1.8).

Results

There were no differences at baseline between controls and patients with incident radiographic knee OA in their clinical covariates (age, gender and BMI, in Table 1) for any of the three definitions (∆JSN, ∆TOS, ∆KL) tested.

Table 2 presents the details of the significant associations of the parsimonious models' features with knee OA initiation for the three definitions. The clustering effects due to the inclusion of both knees were tested using the generalized estimating equations (GEE) [START_REF] Liang | Longitudinal data analysis using generalized linear models[END_REF] . The clusters were defined using the patient ID from the OAI dataset. These results show that considering both knees from each subject helps to better explain the coefficients of the models (JSN, TOS, KL), as parameters considered not significant with a naïve interpretation become significant.

Joint Space Narrowing

The clinical covariates alone were not able to predict the increased JSN. Individual TBT parameters computed in the vertical direction increased the results (AUC > 0.6), yet the prediction was not suitable to discriminate patients with a higher risk of knee OA onset (DOR < 1). TBTp computed in the horizontal direction enabled the prediction (AUC 0.7), however, the prediction was found relevant (DOR > 1) only for the larger scales (m-scale). The combination of all the texture information provided the best prediction (AUC 0.73) with a higher risk for the patients highlighted as at risk (DOR 4.51).

Tibial Osteophytes

The clinical covariates alone were not able to predict the occurrence of TOS. Individual TBT parameters computed in either vertical or horizontal direction increased the results (AUC > 0.6) with relevant potential selection of the patients at risk (DOR > 1) except for the larger scale (m-scale) in vertical orientation. The combination of all the texture information provided the best prediction (AUC 0.71) with a higher risk for the patients highlighted as at risk (DOR 7.10).

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 11

Radiographic knee OA

The clinical covariates alone were not able to predict the occurrence of radiographic knee OA (KL ≥ 2). Individual TBT parameters computed in either vertical or horizontal direction increased the results (AUC > 0.6), yet the prediction was relevant (DOR > 1) only for the horizontal analysis. The combination of all the texture information provided the best prediction (AUC 0.69) with a higher risk for the patients highlighted as at risk (DOR 2.33).

Discussion

In this study, we examined the ability of baseline TBTp measured on computed radiographs in patients from the large cohort of the OAI to predict the incidence of radiographic knee OA over a period of four years. Our findings confirmed that TBTp on baseline radiographs predicted later onset of radiographic knee OA.

Recently, Podsiadlo et al. 19 reported that directional fractal signature analysis was associated Our approach included the whole proximal tibial trabecular bone, mapping the area with 16 different ROIs to cover the global aspect of knee OA as defined by Loeser [START_REF] Loeser | Osteoarthritis: A disease of the joint as an organ[END_REF] . We have already demonstrated that relevant information to detect the risk of OA progression is provided by both lateral and distal subchondral bone regions at the tibia and not only by the medial one [START_REF] Janvier | Subchondral tibial bone texture analysis predicts knee osteoarthritis progression: data from the Osteoarthritis Initiative[END_REF][START_REF] Woloszynski | A signature dissimilarity measure for trabecular bone texture in knee radiographs[END_REF] .

There were no significant differences between the common covariates (age, gender, BMI) at baseline between patients with incident knee OA and controls at 4 years. BMI has been shown to be positively associated with the onset of knee OA [START_REF] Manek | The association of body mass index and osteoarthritis of the knee joint: An examination of genetic and environmental influences[END_REF][START_REF] Jiang | Body mass index and susceptibility to knee osteoarthritis: a systematic review and meta-analysis[END_REF][START_REF] King | Obesity & osteoarthritis[END_REF] . However, BMI seems to be more associated with knee pain [START_REF] Rogers | The association of BMI and knee pain among persons with radiographic knee osteoarthritis: A cross-sectional study[END_REF] than with joint failure [START_REF] Graverand | Progressive increase in body mass index is not associated with a progressive increase in joint space narrowing in obese women with osteoarthritis of the knee[END_REF] . Conversely, we showed that the use of the TBTp could enable the prediction of incident radiographic knee OA.

Details of the TBT analysis

The differentiation between three incidence definitions (∆JSN, ∆TOS and ∆KL) and the dual scale analysis showed that different results are obtained depending on the scenario.

First, it seems that the early changes induced by knee OA are observable along the horizontal direction as reported by Buckland-Wright. These results are in line with previous data shown by Podsiadlo et al. 19 where they found better results with a larger scale in computed radiographs. In our previous study [START_REF] Janvier | Subchondral tibial bone texture analysis predicts knee osteoarthritis progression: data from the Osteoarthritis Initiative[END_REF] we obtained better results when running the analyses in

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

13 the vertical direction for the progression of pre-existing knee OA. This implies that TBT analysis is relevant in the horizontal direction for early diagnosis and then in the vertical one for late follow-up. Extending this observation to the actual topology of the microarchitecture of the trabecular bone implies that the OA process involves changes in the vertical trabeculae (characterized by horizontal analysis) in the early stages of knee OA incidence, while in later stages of knee OA progression these changes may occur mostly on the horizontal trabeculae (characterized with a vertical analysis).

Second, the bifractal aspect of the TBT was arbitrarily chosen and delimited in this study according to empirical observation of the variograms. However, the cut-off at 500 µm corresponds to the breaking point in scale range between the trabecular thickness (100-300 µm) and the trabecular separation (200-2000 µm) [START_REF] Griffith | Bone mass and architecture determination: state of the art[END_REF] . This observation has been reported in several studies focusing on the TBT [START_REF] Kraus | Trabecular morphometry by fractal signature analysis is a novel marker of osteoarthritis progression[END_REF]17,[START_REF] Harrar | Piecewise Whittle estimator for trabecular bone radiograph characterization[END_REF] and might confirm the need to differentiate the process occurring at different scales in bone remodeling. Moreover, the decomposition into partial models for the different definitions showed, as for the orientations, several differences.

The TOS seems to have more links with the global TBT features as three of the four parameters provided relevant results, whereas the higher impact of large scale changes for JSN prediction tends to indicate that the decrease in the joint space width is linked to the organization (m-scale) of the vertical trabeculae.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 14

Strengths and limitations

This study presents several important strengths. The database used included 319 subjects with publicly available radiographs and grading. Both control knees and incident cases had KL 0 at baseline, so our TBT findings are likely representative of early subchondral bone changes that occur before definite radiological OA features are detectable on plain radiographs either in the medial or the lateral compartment. Also, this dataset included multicenter radiographs using different equipment without calibration for TBTp. Neither the variations in the acquisition parameters nor the impact of the soft tissue in the X-ray absorption and diffusion were taken into account. In other words, the descriptors were not adjusted for the equipment or clinical center and the prediction reflected the results of a multicenter cohort without adjustments. Regression models were reduced parsimoniously and cross-validated to avoid highly complex or over-fitted models. Some limitations have also to be considered. The compartments affected with knee OA according to the local TBT were not treated independently and it might be relevant to stratify the 'cases' into several groups (e.g. lateral/medial initiators). Moreover, the femoral subchondral bone texture was not assessed although it might be relevant in this task. In the subject inclusion criteria, digitized films were rejected due to the non-standardized protocol used for digitization. The piecewise fractal aspect of the TBT and a fortiori the cut-off scale have not been completely investigated in this study. However, deeper insight of the

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

15 multifractal characterization of the TBT might be relevant for early OA diagnosis. Finally, the inclusion of the global TBT was performed using a patchwork of ROIs, yet the interactions between the TBTp were not investigated. This requires further work on the characterization of this global TBT with a more parsimonious and meaningful descriptor than a 64-feature vector.

This study highlighted the high selectivity of the actual at-risk patients using the trabecular bone texture. The standard accuracy or AUC values reported in the literature does not reflect details of the sensitivity and specificity of the predictive models. Herein, high values of DOR clearly indicates that TBT analyses can identify a subset of true OA initiators with a very low number of false positives.

Conclusion

Our findings, along with those reported by Podsiadlo et al. 19 , highlight the crucial role of the subchondral bone in the early stage of knee OA. The TBT analysis proposed in this study demonstrated its potential to predict the incidence of radiographic knee OA and highlighted several sub-clusters of texture features (orientations and scales).

The quantification of the TBT of the subchondral bone at the tibia needs further investigation to clearly state the links with the actual trabecular microarchitecture. However, it may contribute to the development of a new and feasible imaging biomarker able to predict disease initiation and identify patients at high risk. This predictive ability would be highly Table 2: Details of the significant associations of texture parameters and covariates with knee OA initiation defined using the joint space narrowing (JSN), tibial osteophytes (TOS)

and Kellgren-Lawrence (KL) grades. These features were extracted from the models reduced with a stepwise algorithm based on the Akaike parsimony criterion.

Table 3: Details of the AUROC and DOR values with 95% confidence interval for radiographic knee OA incidence prediction using the three knee OA initiation definitions based on the joint space narrowing (JSN), tibial osteophytes (TOS) and Kellgren-Lawrence (KL) grades.

Figure 4

 4 Figure 4 presents the ROC curves averaged for the 300 cross-validated runs and Table 2

  with the incidence of radiographic knee OA. The work by Podsiadlo and ours demonstrated the potential interest of selecting a subgroup of patients having baseline characteristic features of TBTp on plain radiographs prone to the development of radiographic knee OA. In their work, 2052 knees were investigated with a KL grade ≤ 1 from 2 centers analyzed individually, focusing on a single pair of ROIs using several fractal descriptors. However, they included OA doubtful subjects (KL grade 1) in their baseline cohort. In the present M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 12 work, only the non-doubtful healthy subjects at baseline were included, as advocated by the OAI 21 .

  E P T E D useful in clinical trials to select patients with the highest risk of initiation or progression of knee OA when new disease-modifying treatments become available. T, Podsiadlo P, Stachowiak GW, Kurzynski M, Lohmander LS, Englund M. Prediction of progression of radiographic knee osteoarthritis using tibial trabecular bone texture. Arthritis Rheum. 2012;64(3):688-695. doi:10.1002/art.

Table 1 :

 1 Characteristics of the 319 patients -344 knees included in this study. ∆JSN: the decrease in the joint space width evaluated as an increase in the OARSI JSN grade over 48 months (in either the medial or lateral side); ∆TOS: the onset of tibial osteophytes evaluated as an increase in the OARSI TOS grade over 48 months (in either the medial or lateral side); ∆KL: the increase in the KL grade with certainty of OA presence after 48 months (KL ≥ 2).

Figure 1 :

 1 Figure 1: Flowchart of the initial OAI cohort pruning (n is the number of subjects and k is

Figure 2 :

 2 Figure 2: Knee radiograph mapped with a customized patchwork of ROIs.

Figure 3 :

 3 Figure 3: TBT variograms computed along the two main loading directions (0° and 90°,

Figure 4 :

 4 Figure 4: ROC curves obtained for the radiographic knee OA incidence prediction using a

*

  the statistical significance of the coefficients was evaluated assuming independent observations (naïve) and accounting for the clustering of knees within a person (robust) using generalized estimation equations.Feature format: "{XY} -{µ;m} | {0°;90°}" {XY} is the ROI identifier in the patchwork (Figure2) X: column from 1 (subchondral) to 4 (deep trabecular) Y: row from 1 (medial) to 7 (lateral) {µ;m} correspond to the scale of observation {0°;90°} correspond to the orientation considered

Area under the receiver-operating characteristic curves (AUROC) Diagnostic odds ratios (DOR) H

  0°/90° correspond to the orientation considered TBTp includes the combination of {H μ|90°; H μ|0°; H m|90°; H m|0°} * combined with the clinical covariates (Age, gender and BMI)

			ACCEPTED MANUSCRIPT
		ΔJSN	ΔTOS	ΔKL
	Cov	0.59 (0.50-0.68)	0.59 (0.51-0.67)	0.57 (0.50-0.64)
	*H μ|90°0 .60 (0.51-0.70)	0.64 (0.55-0.72)	0.60 (0.52-0.67)
	*H m|90°0 .64 (0.55-0.73)	0.61 (0.53-0.69)	0.62 (0.54-0.69)
	*H μ|0°0 .67 (0.59-0.75)	0.66 (0.58-0.73)	0.66 (0.59-0.72)
	*H m|0°0 .69 (0.61-0.77)	0.67 (0.59-0.75)	0.66 (0.59-0.72)
	*TBTp	0.73 (0.66-0.80)	0.71 (0.64-0.79)	0.69 (0.63-0.76)
		ΔJSN	ΔTOS	ΔKL
	Cov	NA	NA	NA
	*H μ|90°0 .01 (0.00-0.19)	3.46 (3.13-3.83)	0.45 (0.36-0.56)
	*H m|90°0 .68 (0.51-0.90)	0.06 (0.00-0.60)	0.06 (0.04-0.09)
	*H μ|0°0 .85 (0.68-1.04)	5.71 (5.21-6.21)	1.78 (1.66-1.91)
	*H m|0°5 .17 (4.72-5.68)	5.84 (5.46-6.29)	3.49 (3.32-3.67)
	*TBTp	4.51 (4.27-4.77)	7.10 (6.67-7.43)	2.33 (2.24-2.43)

is a fractal parameter linked to the complexity µ/m correspond to the scale of observation NA: not applicable, values could not be computed due to high misclassification rates

Acknowledgments

The authors would like to thank the OAI study participants and clinical staff as well as the coordinating center at UCSF. The authors wish to thank the ANR for financial support under the project ANR-12-TECS-0016-01. We acknowledge Mrs Villequenault for typing the manuscript.

OAI disclaimer

The OAI is a public-private partnership comprising five contracts (N01-AR- Private sector funding for the OAI is managed by the Foundation for the National Institutes of Health. This manuscript was prepared using an OAI public use data set and does not necessarily reflect the opinions or views of the OAI investigators, the NIH, or the private funding partners.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Author Contributions

All authors contributed to the conception and design of the study, analysis and interpretation of data, drafting the article and revising it critically for important intellectual content, and approved the final version submitted. Dr Lespessailles takes responsibility for the integrity of the work as a whole, from inception to finished article.

Role of the funding sources

Thomas Janvier (PhD student) was funded by ANR-12-TECS-0016-01.

Funding sources had no role in the study design, collection, analysis and interpretation of the data or the decision to submit the manuscript for publication.

Conflict of interest

All authors state that they have no conflicts of interest.