Michael Morin

Rallou Thomopoulos

Irène Abi-Zeid

Maxime Léger

François Grondin

Martin Pleau

Explaining the Results of an Optimization-Based Decision Support System -A Machine Learning Approach

In this paper, we present work conducted in order to explain the results of a commercial software used for real-time decision support for the flow management of a combined wastewater network. This tool is deployed in many major cities and is used on a daily basis. We apply decision trees to build rules for classifying and interpreting the solutions of the optimization model. Our main goal is to build a classifier that would help a user understand why a proposed solution is good and why other solutions are worse. We demonstrate the feasibility of the approach to our industrial application by generating a large dataset of feasible solutions and classifying them as satisfactory or unsatisfactory based on whether the objective function is a certain percentage higher than the optimal (minimum) objective. We evaluate the performance of the learned classifier on unseen examples. Our results show that our approach is very promising according to reactions from analysts and potential users.

Introduction

Decision support systems (DSS) and in particular optimization systems are more and more common in engineering, business, and management environments. Although these systems may be very sophisticated and powerful in terms of solving complex mathematical models and providing optimal recommendations, more often than not, the users of these systems are skeptical of the results presented. This leads to a resistance and a mistrust with regards to these technologies. More than ever, it is becoming important to explain the solutions of the optimization systems in terms that users and operators can understand, in order to be comfortable in following or challenging the recommendations of these systems. In this paper, we present work conducted in order to explain the results of a commercial software used for real-time decision support for the flow management of a combined wastewater network. This tool is deployed in many major cities and is used on a daily basis. Over the years, our industrial partner has identified the need for explanations of the optimization results which prompted us to define a collaborative project.

Du et al. [START_REF] Du | An explanation oriented dialogue approach and its application to wicked planning problems[END_REF] distinguish between "backward" and "forward" explanations. Backward explanations are past-oriented. They seek reasons in the past (actions, events) that led to the observed situation. Forward explanations are future-oriented. They foresee the future consequences of a possible decision. In the present paper, we propose a forward approach which provides explanations for the recommendations of an optimization-based DSS. We apply a machine learning technique to build rules for classifying and interpreting the solutions of the optimization model, namely decision trees. Our main goal is to build a classifier that would help a user (not a researcher or a practitioner from a field related to optimization) to understand why a proposed solution is good and why other solutions are worse.

The paper is structured as follows: Section 2 presents a brief literature review, Section 3 includes the problem background and methodology, Section 4 contains a description of the experiment and the results. We conclude in Section 5.

Literature Review

Decision support systems have been around for several decades now. Interestingly, since the 90s, several studies have pointed out that end users have a low acceptance level of the solutions computed by DSSs [START_REF] Davis | Determinants of decision rule use in a production planning task[END_REF][START_REF] Kleinmuntz | Why we still use our heads instead of formulas: Toward an integrative approach[END_REF][START_REF] Lehtola | Requirements prioritization challenges in practice[END_REF][START_REF] Lim | Judgmental forecasting with time series and causal information[END_REF]. These studies have explored, respectively, the importance of feedback about a decision made, the role of cognitive biases in human decision errors, the complexity of the requirements prioritization which implies an iterative decision-making process, and the importance of contextual and causal information in forecasting.

However, during the 2000s, studies about the effects of providing explanations on DSS advice acceptance were still extremely rare [START_REF] Lawrence | Improving forecast utilization by providing explanations[END_REF] and they have only recently begun to emerge as an important consideration [START_REF] Du | Two machine-learning techniques for mining solutions of the ReleasePlanner™ decision support system[END_REF]. In their paper, Du et al. [START_REF] Du | Two machine-learning techniques for mining solutions of the ReleasePlanner™ decision support system[END_REF] use two machine-learning techniques, namely rough sets and dependency networks, to mine DSS solutions. The mining results are then provided as explanations of the DSS solutions. Both techniques provide different, but not disjoint, results, which are used by the authors to classify the mining results as strong or weak explanations.

In [START_REF] Lawrence | Improving forecast utilization by providing explanations[END_REF], two types of explanations are distinguished: technical ones, tracing the process used to make a forecast, and managerial ones, providing a justification based on the meaning of the forecast. Both types of explanations were estimated to increase decision-makers confidence. This diagnostic is confirmed by Gonul et al. [START_REF] Gönül | The effects of structural characteristics of explanations on use of a DSS[END_REF]. Based on their review and conclusions, we have built the following three tables which summarize some conclusions about the categories, the roles and the added value of explanations. These are presented respectively in Tables 1, 2 and 3.

Background and Methodology

The commercial DSS that we are interested in is a software for real time control and management of a combined wastewater network, provided by our industrial partner. It includes a simulator and an optimizer. The network is made-up of controlled and uncontrolled sections. Based on weather forecasts over a given time period, the simulator simulates the flow in the uncontrolled parts of the network that provide the inputs to the controlled parts of the network. The optimizer's role is to provide the complete set of control point configurations (decision variables values) that minimizes the overflow cost (objective function) for the considered period. It establishes, based on water flows in the network, a set of instructions that are sent to local stations to regulate flows using sluice gates and pumps. Constraints describe the behavior of the network and the cost function is defined to minimize overflow volumes, storage facility dewatering time, to ensure a balanced hydraulic load distribution throughout the network and to minimize the associated management costs [START_REF] Pleau | Global optimal real-time control of the Quebec urban drainage system[END_REF]. The provided solution,

Effects on users Consequences Result

Better understanding Learning

Perform better decision making (accuracy and speed) Accept the system as logical and thus grounded and useful Follow the logical reasoning used by the system Positive perception of the system (ease of use, usefulness, satisfaction, confidence) although optimal with respect to the objective function in the model used by the optimizer, can be modified by the end-user.

A solution fully describes, in terms of the values of the decision variables, the state of the controlled section of the network. Decision variables determine, in particular (but not exclusively), the positions of all the control gates regulating the wastewater system flow management. These can be closed or open to a certain extent.

Methodology

Given an optimal solution to an optimization problem provided by solving a mathematical program, in this case a minimization problem, we wish to enable an end-user to understand which characteristics of a solution render it more or less suitable (satisfactory) from the perspective of the decision s/he has to make. From a single objective optimization perspective, a solution that is optimal minimizes a predefined and accepted objective function while being feasible, i.e., while respecting a set of predefined constraints. We employ the term satisfactory to describe a feasible solution that is a suitable decision in the context where a mathematical program is used as part of a DSS. An optimal solution is of course satisfactory. However, some satisfactory solutions might be suboptimal.

Suppose, for instance, that a small increment δ in the objective function value is tolerable. Any feasible solution x with an objective value f (x) that is not worse than f (x * ) + δ , i.e., f (x) ≤ f (x * ) + δ , is said to be satisfactory. Although x is not optimal, it is a suitable decision with respect to that criterion.

We therefore need a methodology to characterize the satisfactory solution set, using simple rules, in order to allow a DSS user to understand the following:

1. why a solution (including the optimal one) is satisfactory; and 2. why a given change to an optimal solution might produce a solution that is worse, or potentially catastrophic in practice.

If a user is not satisfied with the optimal solution initially provided by the DSS for a number of possible reasons, including the fact that it does not intuitively make sense to her/him, the rules (explanations) would allow her/him to understand which decision variables can be modified (and how) before the DSS can re-optimize and return a new satisfactory solution.

We chose to use decision trees to produce (learn) such a simple set of rules to characterize the set of satisfactory solutions. This is similar to the approaches in [START_REF] Du | Two machine-learning techniques for mining solutions of the ReleasePlanner™ decision support system[END_REF], albeit in a different context and on bigger-size data sets.

Decision trees and the classification of satisfactory solutions

Decision trees are a widely used machine learning technique. In the next subsections, we briefly present machine learning and decision trees as we applied them to the problem of characterizing the set of satisfactory solutions.

Machine learning

(Supervised) machine learning is the art of learning by examples [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF]. A learning algorithm is provided with a training set S of examples. Each example in set S is a pair (x, y) where x is an input and y is an output. The output can be a class, a real number, or it can have a structure [START_REF] Bakȋr | Predicting Structured Data[END_REF]. When x ∈ R n and y ∈ {True, False}, we have a binary classification problem. We call the entries in x the features and y the class.

A solution to the binary classification problem is a function h, called a classifier, that takes an input x ∈ R n and that outputs y ∈ {True, False}. A learning algorithm builds a classifier h using the set S, i.e., known pairs (x, y) of features and class (training phase). Once built, the classifier can be used to compute an output class y given an unseen vector of input features x (testing phase). In this paper, we formulate the problem of recognizing satisfactory solutions from unsatisfactory ones as a binary classification problem.

Decision trees

A decision tree [START_REF] Breiman | Classification and Regression Trees[END_REF] is a classifier in the form a tree structure. Classifiers can be either white box or black box models. We chose to use decision trees since a decision tree is a white box and interpretable model made of simple decision rules. It takes x as an input, and it outputs a class y. The nodes of the tree are either decision nodes or leaves. A decision node corresponds to a branching rule based on the value of a single input feature. A leaf node is assigned a prediction class.

In Other stopping criteria to manage the height of the tree can be used. In that case, not all the examples in a partition that reach a leaf will have the same class and the value of the leaf can be set to the majority class.

In the testing phase, an unseen example x, i.e., an example for which y is unknown, follows a path from the root to a leaf. The example is tested against each decision rule. If it respects the rule, it branches left, otherwise it branches right. The predicted class is that of the reached leaf.

Learning to recognize satisfactory solutions

As mentioned previously, each node of a decision tree, except the leaves, is a branching rule. A path from the root to a leaf can be seen as a rule describing the satisfactory solution set. We first learn a decision tree (training phase). Later on, this decision tree is used to determine whether a solution is satisfactory or not (test phase). Furthermore, since a decision tree is interpretable, it provides a set of rules characterizing satisfactory solutions. We next describe the process used to generate the set S of training examples that can be used for learning.

Assume that an optimal solution has been provided by the DSS and that a satisfactory criterion has been determined. We generate set S by simulating a sufficiently large number of modifications of the optimal solution. Simulating a modification consists in randomly fixing a decision variable's value. After each simulation, we evaluate the value of the solution according to the chosen criterion that discriminates between satisfactory and unsatisfactory ones. We store in S, the value of all decision variables along with the class of the newly generated solution: satisfactory (true) or not (false). Finally, we learn one or more decision trees on S.

Experiment

In order to demonstrate the feasibility of the approach to our industrial application, we generated a large dataset of feasible solutions and classified them as satisfactory or not based on whether the objective function is a certain percentage higher than the optimal objective. We then used a small subset of these solutions to learn a classifier. The small subset used for learning is the training set (S). We evaluated the performance of the learned classifier on the remaining unseen examples. We provide the details on the data we generated for the experiment in the next subsection.

Data preparation

In order to learn the characteristics of the satisfactory solution set, we first obtained an optimal solution to a problem instance using the optimizer of the DSS developed by our partner. The weather forecast, the simulation results, and the provided optimal solution belong to a real use-case scenario.

Starting from the optimal solution, we generated 10,000 feasible solutions by simulating a single modification of one of the decision variables by the end-user. Each solution of the dataset used for the experiment is generated as follows:

• one of the decision variables is chosen and modified randomly;

• the objective value of the solution resulting from the perturbation of that single control point is reevaluated.

The resulting dataset contains the values of all the variables of the optimization model for all the generated solutions along with their objective value. We then filter the dataset to keep only the variables that can have an explanatory significance for the end-user. The identification of these variables was performed by the domain experts. We excluded all the variables that are used as components of the objective function. We then assigned a class to each solution. For the purpose of the experiment, we assume that a solution is nonsatisfactory if and only if its objective value exceeds by 5% (or more) the optimal value. In the next subsection, we present, the results of our methodology when used jointly with decision trees.

Decision trees and their predictive accuracy

All decision tree experiments are conducted using scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF], a machine learning library for the Python programming language.

Figure 1 shows the learning curve we obtain on our industrial instance. Most learning algorithms for decision trees are greedy algorithms. Scikit-learn's learning algorithm for decision trees makes no The decision rules learned by the decision trees were validated by experts from the application field. The tendency of the learning algorithm was to choose high impact control points (decision variables) in its first few decision rules.

Conclusion

In this study, we have proposed an approach using decision trees to explain the results of a commercial software used for real-time decision support for combined wastewater system flow management. Our solution allows analysts to identify the main criteria that influence the DSS outcome, and helps users in understanding the rationale behind the DSS's recommendations. The rationale supporting this approach consists in classifying the space of solutions according to some criterion that characterizes the performance of the solutions proposed by the DSS.

We plan to explore several directions in the next stages of this work. First, we will identify "optimal" sizes of the training set and tree depths, allowing for both a good learning quality and a reasonable computing time. Moreover, the robustness of the method will be tested based on the following factors: on the one hand, the number of variables taken into account, in particular for increased numbers of predictive variables; on the other hand, the proportion of positive examples composing the training set. Finally, we plan to experiment our approach on other kinds of learning algorithms, also able to compute interpretable classifiers, with the long-term objective of providing users with arguments expressing the learnt rules.

  the training phase, a decision tree is built in a greedy fashion by the learning algorithm. The examples of the training set S, i.e., the known examples, start at the root of the tree. A first branching rule is applied at the root on the training examples. The learning algorithm chooses the branching rule that best partitions the training examples according to a heuristic. All examples x that respect the branching rules take the left path. All other examples take the right path. At the next level, a new decision node is created for each partition and two branching rules are selected according to the partitioning heuristic. Branching stops when all the examples of a partition belong to the same class.

Figure 1 .

 1 Figure 1. Learning curve: Mean accuracy of the learned decision trees (solid curve) and of the simple classifiers (dashed curve) on the test set against various training set sizes with standard deviation in blue and green shading and upper and lower bounds in light blue and light green

Table 1 .

 1 Categories of explanations

	Categories	Subcategories	Remarks
	Depending on the content of	Explanations on the how: system
	explanations	reasoning mechanism (information
		used, rules, steps)
		Explanations on the why:	The most efficient for the
		justification of a decision	acceptance of the result
		highlighting the reasons behind the	provided
		decision
		Strategic explanations: overview of	Less used since their
		the whole issue and the stakes	obtention is difficult
		tackled by the system
	Depending on the way of	On user demand
	providing explanations	Automated
		Intelligent: when a need for	Complex mechanism
		explanations is detected by the
		system
	Depending on the format in	Text-based: predefined sentences	Natural language-style
	which explanations are	or combinations following some	sentences improve
	presented	rules	understanding and acceptance
			of the system
		Multimedia: graphs, pictures,	Costly implementation but
		animations	persuasive and efficient to
			obtain users' confidence

Table 2 .

 2 Roles of explanations

	Objectives	Benefits for users
	Explain a detected anomaly	Understand a recommendation
		Verify that it satisfies the expectations
		Solve a contradiction between the users and the system
	Provide complementary information	Participate in solving a problem (short-term)

Facilitate learning

Perform better in future tasks (long-term)

Table 3 .

 3 Added value of explanations

e-mail: Michael.Morin.3@ulaval.ca e-mail: Irene.Abi-Zeid@osd.ulaval.ca


Grégoire Lenoble
email: gregoire.lenoble@siemens.com

Marie Olivier

David Donjat

Anton Steblinkin

Contribution de l'analyse thermique MBSE et 3D à la conception du moteur de l'actionneur de stabilisateur d'un avion régional

Résumé -Aircraft hydraulic systems traditionally provide power for flight surface actuation, but are heavy, complex to maintain in service, and consume significant amounts of energy. Electrical based actuators help tackle these issues, but come with other challenges, such as thermal integration. Indeed, the close vicinity of several heat sources (hydraulic pump and bloc, electrical motor, control unit) in a cavity that is only passively cooled makes it difficult to guarantee the electrical motor and Electronic Control Unit (ECU) safe operation over the entire aircraft mission from a thermal point of view. In addition, systems thermal behavior can be critical when used inside composite structures that can lose durability under high temperatures. Collaborative European-Russian project RESEARCH (w3.onera.fr/RESEARCH/) has enabled the thermal analysis of a regional aircraft stabilizer bay that houses two elevator actuators, so as to support and optimize actuator motor design. This analysis largely relied on Model Based Systems Engineering (MBSE) given the early actuator design phase and the need for mission-wide analysis, through a number of varying scenarios. However the additional insight of three-dimensional analysis has also been included with simulations performed at thermally critical flight conditions. The methodology, models, outputs and implications for motor design gained from this study are described in this paper.

Nomenclature

Introduction

The reduction of fuel burn has become a priority for most aircraft operators due to environmental and cost implications. Such interest has encouraged the aerospace industry to look at ways of minimizing fuel consumption, concluding that the approach towards a more electrical aircraft, i.e. using electricity to operate the consumers, could provide some benefits. This approach tends to remove hydraulic actuators and the hydraulic systems that feed them, as they require a significant amount of maintenance effort, are heavier and consume significant amounts of energy. Additionally, this kind of electrical components allow for a higher integration of sensors for the monitoring of actuators performance and structural integrity (from a safety point of view); or achieving an improved degree of "intelligence", improving the actuation system integration within the controlling system or aircraft. Despite the potential benefits indicated above, such as weight reduction, increased safety and reliability, lower consumption, reduced maintenance costs, etc… electrical actuators are not fully embraced by the industry due to some of the limitations they present in their current development stage and which prevent them from being used massively or in critical areas, such as the aircraft primary flight surfaces. These limitations are: high actuator weight, safety concerns due to the jamming probability of electromechanical actuators (EMAs) and Electro Magnetic Compatibility (EMC) issues. Systems thermal behavior and integration is one of the most challenging among these issues. For this reason, and the RESEARCH program being a wide objective oriented cooperative project, the decisions was taken to also cover the task of managing the thermal integration of the actuation system under development by the consortium partners.

It's well known that the reliability of EMA or electro-hydrostatic actuators (EHA) systems for flight controls is dependent on their operating temperature. The A380 flight control EHA is a well-done example [START_REF] Van Den Bossche | The A380 flight control electrohydrostatic actuators, achievements and lessons learnt[END_REF] as a backup actuator. But EHA and EMA constantly working in active regime need to be improved from a thermal management viewpoint, as the elimination of the classical hydraulic network imposes the need to dissipate heat generated by the actuation system (motor and ECU, and also harness) differently. The management of this thermal environment needs to develop both new cooling technologies (see [START_REF] Lawson | Thermal Management of Electromechanical Actuation On An All-Electric Aircraft[END_REF]) and design optimization methods. This last improvement target appears as a key point when considering thermal aspects earlier in the design process.

Cavity configuration

Based on the results of weight and safety analyses of the different possible FCAS architectures performed by the RESEARCH consortium in the first project phase, the elevator control system of a regional passenger aircraft (100-seater) was chosen as the target application to be developed. The control system should consist of an EMA (under CESA's development, Spain) and EHA (under Moscow Aviation Institute (MAI) development, Russia) operating in parallel (active/active configuration) (see figure 1 with EHA configuration). The overall idea was to optimize the actuators design in terms of acceptable heat losses and the heat dissipation taking into account the actuator real duty cycle and cavity geometry at the early design stage.

The flight mission

The methodology of the actuation system thermal behavior analysis described further implies the wide use of simulation, which as a first step was represented by the simulation of the cavity system thermal behavior modeled in LMS Amesim during the typical flight of a regional jet. The flight scenario is based on the real flight record and contains the data arrays of elevator actuator displacement and external forces under hinge moments as well as other flight parameters (altitude, true air speed etc.), which thus fully determines the actuation system heat losses and heat dissipation processes. Also in order to obtain the duty cycle of the actuators during some highly loaded flight scenario, a «virtual» aircraft flight-record was also developed. It contains the influence of predefined turbulence according to Figure 2. This scenario also includes some special maneuvers during flight (e.g. s-maneuver just before landing) to put the actuation system into the most critical regime.

Thermal analysis 2.1. MBSE cavity model

Cavity boundary conditions

The elevator actuator operates inside a given thermal context, which will interact with the actuator bay, and even influence its thermal behavior. The following interfaces to the actuator are considered: outside air, Sun radiation, cavity structure. Outside airflow provides a large cooling effect to the cavity outside shell. This is captured both on the top and bottom panels, assuming airflow over a horizontal flat plate configuration (see Figure 3). Depending on the flight regime and air properties, forced and free convections are calculated with laminar or turbulent equations [START_REF] Incropera | Fundamentals of Heat and Mass Transfer[END_REF]. In order to ensure a smooth mathematical model over the entire range of convection regimes (purely free, mixed, purely forced), a unique Nusselt number combines the free and forced ones, so that a unique convection Heat Transfer Coefficient (HTC) can be calculated as per [START_REF] Van Den Bossche | The A380 flight control electrohydrostatic actuators, achievements and lessons learnt[END_REF].

𝐻𝑇𝐶 = (𝑁𝑢 * 𝜆)/𝐿 (1)
In this outside convection case, the characteristic dimension L is the total stabilizer outside area (top and bottom) divided by its perimeter. Note that the air viscosity, thermal conductivity, specific heat, density and speed are dynamically updated during the aircraft mission simulation.

As shown with Figure 3, Sun radiation onto the stabilizer cavity is accounted for with (2), where θ is the angle between the cavity surface (normal vector) and the Sun direction (set to 0° so as to consider maximum Sun radiation), and abscoeff stands for the surface absorption factor.

𝜑 𝑆𝑜𝑙𝑎𝑟 𝑝𝑎𝑛𝑒𝑙 = 𝜑 𝑠𝑜𝑙𝑎𝑟 × cos 𝜃 × 𝑎𝑟𝑒𝑎 × 𝑎𝑏𝑠𝑐𝑜𝑒𝑓𝑓

This Sun radiation is complemented with the radiative exchange the surface has with outside ambient air, as per (3), and using ε gw as the outside air and cavity panel emission factor.

𝜑 𝑟𝑎𝑑 = 𝜎 ⋅ 𝜀 𝑔𝑤 ⋅ 𝑎𝑟𝑒𝑎 ⋅ (𝑇 𝑠𝑢𝑟𝑓𝑎𝑐𝑒 4 -𝑇 𝑎𝑖𝑟 4 ) (3)
Note that equation ( 3) is a particular case of the general radiation equation that assumes an infinite "radiating air surface" and thus a view factor of 1.

As shown with Figure 4, ground radiation is accounted for in a very similar way, only with no exposure angle factor, a different ε and the ability for it to be turned off when the aircraft has left ground. The thermal behavior of the cavity that houses the stabilizer actuator has been simplified to that of a few panels: a vertical, and two horizontal ones: top and bottom. Top and bottom panels are the aggregate of an external layer (exposed to the Sun / ground and outside air), and an internal layer (in contact with internal air), separated by an insulation composite material, as shown with Figure 3 and Figure 4.

Finally, thermal conduction occurs from the top to the bottom panels, through the vertical one. This is captured with two instances of a "thermal conduction" component, that derive the heat flux between two solid material pieces 1 and 2; depending on is their contact surface, and contact thermal resistance. It predicts T1 and T2 temperatures, both taken at a distance d1 and d2 inside each material piece.

Cavity internal heat loads

The end objective of this analysis being the thermal integration of all heat sources and exchanges inside this cavity, each heat source has then been modeled individually. The heat transfer modes involved in all of the following heat loads are all based on standard heat transfer equations: free and forced convection over a flat plate, radiation in an enclosure, and linear conduction between two pieces of different material.

In this article, we are focusing on EHA systems which are modeled with a focus on their thermal interaction with the cavity. The EHA heat flux generated as a function of its duty cycle is obtained from the EHA Simulink-model developed by MAI, and split over the EHA parts: body, motor, pump, piston and reverse valve. The EHA parts are in obviously contact, and bolted to the cavity bottom (lower panel). They dissipate EHA heat losses through conduction (to themselves, and the lower panel), forced oil convection (driven by duty cycle oil consumption), as well as air convection (forced and free) and radiation, as will be described in the next section.

The Electronic Control Unit (ECU) adds an additional third heat load into the cavity, through its thermal losses. This is captured with a single thermal capacitance. Its heat load is calculated based on the power it provides to the EHA motor, and its own efficiency η. The ECU is in thermal contact with the lower panel, and is cooled by air convection and radiation.

Cavity thermal integration

Thermal integration of the components described in the above sections is mostly piloted by the cavity internal convection and radiation. First, internal convection thermally connects all cavity elements: EHA / EHA, ECU and cavity panels. As seen with Figures 4 and5, convection is modeled over "horizontal" and "vertical" faces, so as to capture all heat exchanges, distinguishing between the EMA cylindrical body and its flat end faces for example. Natural and forced convections are calculated following equations given in [START_REF] Incropera | Fundamentals of Heat and Mass Transfer[END_REF] with the additional subtlety that convection is piloted from the airspeed point of view. LMS Amesim being a Systems Engineering tool, 2D or 3D airspeed maps cannot be predicted. Rather, this air speed information is extracted from ONERA predictions, and condensed into speed profiles that are driven by the ingested airflow, or artificially set to a minimum of 0.2 m/s in the worst case assumption of a sealed cavity, as recommended in [START_REF] Saury | Mesures de vitesse, de température et de flux thermique en convection naturelle en espace confiné[END_REF]. Internal radiation completes the list of cavity internal heat exchanges, and has been allowed between all bodies (EHA, ECU) and both cavity horizontal panels. Radiations between body and from a panel to the other have been neglected given temperature levels and relative influence with regards to convection or conduction. Finally, all heat loads, heat storage, and heat exchanges are resolved in the Cavity submodel shown in Figure 6. This submodel computes the cavity temperature and Pressure based on mass and energy conservation, respectively owing from airflows (ingested and leaked), and the heat exchanges detailed above.

CFD cavity model

CFD model is based on a 3D discretization of the stabilizer cavity as shown on Figure 1. Calculations are performed with FLUENT14's unsteady incompressible solver. According to our estimation of flow conditions, laminar approach is selected. For thermal problem, both free convective and radiative behaviours are resolved. For radiative approach, S2S method is preferred for 3D enclosed configuration. Boundary conditions are extracted from LMS Amesim software simulation in term of heat loads (i.e. heat power released by electric items) and external airflow (i.e. SAT (or Static Air Temperature), altitude and Mach number).

First results

MBSE model results

The evolution of the temperature of the active EHA over the whole flight mission is shown with Figure 7. As expected, all temperatures clearly increase during both ground phases as a result of medium to high SAT and low aircraft speeds. On the opposite, low outside temperatures and high aircraft velocities during cruise directly help reduce EHA temperatures. The steeper temperature rise during approach and landing than during takeoff is a direct consequence of the quick decrease of external heat convection following landing. The temperature profiles observed for the EHA pump and motor are both higher than that of other EHA parts, but yet exhibit different time variations. This is explained by the fact that they both reject the same thermal power (approx. 30W), but have different weights (as summarized with Table 1), which results in different heat per unit mass values. The watchful reader may observe that the flight mission shown with Figure 2 has been extended by 160mins in the taxi-in phase. Although highly unrealistic, this was done so as to verify long-term temperature stabilizations, and made possible by the very low CPU run times of the LMS Amesim software model (100 CPU s for 16800s simulated). Despite this mission extension, all temperatures keep rising after landing and taxi-in phases, even with moderate heat release levels (6W for the EHA, and 20W for the ECU). This can be explained with the conservative assumption of a sealed cavity, combined with a low wind speed assumption (1 m/s) that do not provide sufficient cooling to the cavity and its actuators. This behavior does not hold when cavity ventilation is introduced, as described in section "EMA design support". The results shown above have confirmed that, as expected, the flight condition that is most critical and demanding from a thermal point of view is the landing / taxiing phase. Therefore, it has been decided to further analyze the cavity thermal behavior straight after landing (t=95min) with the CFD model described previously. The end of cruise condition (t=76,7min) has also been analyzed, given that it is the operating point where actuators are designed to operate for the longest part of their life.

CFD model results

The temperatures of the cavity bodies obtained from CFD analysis in an EHA-EHA configuration, on a HOT day, in steady state conditions, are given below as body averages to allow for a comparison with MBSE results (extracted from Figure 7). End of cruise conditions (t=76,7min) are given and compared with MBSE results in Table 1. Interpretation of the MBSE and CFD results provided in Table 1 is of great interest, providing average differences around 6 to 7% for the values that should be compared, 5% as a minimum discrepancy, and 11% as a maximum one. The main contributors to these discrepancies are believed to be differences in local HTC values, and radiation view factors. In this result post-processing summary percentages, the CFD cavity temperature at End of Cruise is not included, as it clearly highlights a limitation to the current model. This is believed to be linked to the inaccurate accounting of thermal conduction along the cavity walls, currently modelled with the "shell-conduction" approach that uses an artificial meshing based on the real panel thicknesses. Improvements to this approach are currently being investigated. Also left out of the post-processing summary percentages is the cavity upper panel external temperature since the end of cruise difference is 1degC only, but also since solar radiation is not yet accounted for in the CFD model. This will be added in future versions of the model, given its importance when compared with other heat fluxes into the cavity roof, as understood from the MBSE model (250 to 350W, compared to 200W for external convection). Keeping in mind the added value of the CFD model for an accurate prediction of radiation view factors, an additional run has be performed with the CFD model, increasing all emissivity values to ε=0.9 (starting from 0.5 for the glass fiber and aluminum, and 0.8 for the carbon fiber). This additional CFD run provides further understanding of the cavity thermal behavior, and especially of the importance of internal radiation, which increases all internal temperatures by 2-3 °C. Just like the cavity air temperature discrepancy, this variation as a result of additional radiation can be explained with the "shell conduction" model that needs improving. Indeed, in the previous (ε=0.5) run, inaccurate thermal conduction between the EHA constituents prevented realistic spread of the heat dissipated by the motor and pump. The rise in radiation (ε=0.9) reduces this isolation of heat in one EHA constituent only (mostly the pump), and allows for heat propagation to the rest of the cavity and hardware. Apart from conduction and solar radiation to be refined in the CFD model, this results comparison still validates the discretization assumptions made in the MBSE analysis, with the split of EHA actuators into a few bodies only. This is of great importance, as it gives confidence in further use of the MBSE model, whose flexibility and fast simulation capability are instrumental at this early stage of actuator design. Finally, further understanding of the flow physics at stake in the cavity is provided with Figure 9, which shows the temperaturecolored profile of some cavity air flow.This illustration is key to the understanding of the natural convection happening within the cavity. A specific post-treatment will be designed during further work to post-process this information for the update (improvement) of HTC estimations used in LMS Amesim software model, basing them on local air velocities as influenced by the cavity 3D geometry rather than lumped estimates.

Temperature (°C

Conclusions

The proposed methodology obviously gives the opportunity to take into account such a critical aspect for More Electrical Aircraft (MOA) applications as system thermal behavior at the early design stage. As a result, it leads to the reduction of the systems development time and what's probably more important the output product is able to represent the more effective solution in terms of mass and dimensions leading to the reduction of the aircraft lifecycle costs. In addition the methodology used gives the ability to all the project partners to fit the outputs of their work all together on the common LMS Amesim language while focusing on their own area of expertise.

It also should be noted that several further research activities planned to become the logical continuation of the work described in the paper as the analysis of actuation system thermal behavior during highly loaded / emergency flight scenarios with turbulence and different failures and also the experimental validation of the thermal model developed in LMS Amesim and 3D calculations from ONERA. The validation is planned to be done in TsAGI with the use of the full-scale composite stabilizer and the real actuators installed inside the cavity.

  emission factor Thermal conductivity (W/m/K)

Figure 1 :

 1 Figure 1: Stabilizer cavity with actuating devices: EHA/EHA/ECU (ONERA)

Figure 2 :

 2 Figure 2 : «Virtual» highly loaded flight scenario

Figure 3 :

 3 Figure 3 : Cavity top panel

Figure 4 :

 4 Figure 4 : Cavity bottom panel

Figure 5 .

 5 Figure 5. EHA model (left) and ECU model (right)

Figure 6 :

 6 Figure 6 : Cavity model focus

Figure 7 .

 7 Figure 7. Time profiles of the active EHA lumped temperatures on HOT days, for a sealed cavity (no air ingestion)

Figure 9 .

 9 Figure 9. Cavity airflow profile colored by air temperature, HOT day, sealed cavity, taxi-in conditions

Table 1 .

 1 Time profiles of the active EHA lumped temperatures on HOT days, for a sealed cavity (no air ingestion)

	)	CFD MBSE	∆	CFD MBSE	∆
	Mean EHA	51.9	54.3	5%	52	54.9	6%
	ECU	47	52.4	11%	63.2	55.8	12%
	Cavity Upper panel -internal Upper panel -external	15.1 ---3	45 35.6 -4	198% --33%	53.2 --50	53.6 53.2 50	1% --0%

Table 2 .

 2 CFD averaged temperatures, steady state at landing conditions, sealed cavity, HOT day, ε=0.9


Damien Lam

Cyndélia Cladière

Katja Guillaume

Susanne Wassmann

Bolte

France Lam

Damien Cladière

Cyndélia Guillaume

Katja Wassmann

Susanne Bolte
email: susanne.bolte@upmc.fr

Super-resolution

Keywords: Confocal microscopy, deconvolution, refractive index matching, 3D-SIM-like resolution improvement, meiotic spindle

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

During the past ten years one of the most fundamental challenges in 3D imaging was to break the diffraction barrier of light in optical microscopy to attain super-resolution images [START_REF] Schermelleh | A guide to super-resolution fluorescence microscopy[END_REF]. These new super-resolution techniques may be based on the precise localization of fluorochromes, such as PALM/STORM, which needs high power excitation and numerous acquisitions to collect enough photons. On the other hand, non-linear fluorophore responses such as Stimulated Emission Depletion (STED) help increasing the resolution. STED is based on the selective deactivation of fluorophores by a high power depletion laser, minimising the area of illumination at the focal point, and thus enhancing the achievable resolution for a given system. Sample preparation in PALM/STORM and STED-microscopy demands for dyes with specific characteristics and fluorochrome environment, especially if multicolour labelling is needed. Another major drawback is that 3Dinformation is difficult to achieve and the depth penetration of these techniques is limited to a few µm for PALM/STORM and below 20 µm for STED.

3D-Structured Illumination Microscopy (3D-SIM), a technique based on illumination patterning has

recently emerged as a practicable approach to circumvent the diffraction limit of light, avoiding complex sample preparation protocols. 3D-SIM requires the acquisition of at least 15 images per plane and the image contrast has to be good enough to extract high frequencies. 3D-SIM allows obtaining a spatial resolution of 100-150 nm laterally and 250-350 nm axially, thus improving resolution by a factor of approximately two compared to confocal microscopy [START_REF] Schermelleh | A guide to super-resolution fluorescence microscopy[END_REF][START_REF] Gustafsson | Surpassing the lateral resolution limit by a factor of two using structured illumination microscopy[END_REF][START_REF] Fischer | Microscopy in 3D: a biologist's toolbox[END_REF]. This technique has enabled high resolution imaging of the cellular components such as the microtubule (MT) network in cells and organisms [START_REF] Gustafsson | Three-dimensional resolution doubling in wide-field fluorescence microscopy by structured illumination[END_REF][START_REF] Winter | Two-photon instant structured illumination microscopy improves the depth penetration of super-resolution imaging in thick scattering samples[END_REF]. However, standard 3D-SIM is limited to an imaging depth of 10-20 µm [6] and therefore difficult to achieve in thicker samples. The reasons for this are the inherent out-of-focus background of this wide-field approach. Furthermore, bulky samples do not allow obtaining a good fringe contrast necessary for the reconstruction of the 3D-SIM-image. Using two-photon SIM, the depth penetration limit has been broken enabling super-resolution imaging in depth beyond 45 µm [START_REF] Winter | Two-photon instant structured illumination microscopy improves the depth penetration of super-resolution imaging in thick scattering samples[END_REF][START_REF] York | Resolution doubling in live, multicellular organisms via multifocal structured illumination microscopy[END_REF]. Yet, 3D-SIM and 2-photon SIM are methods that call for specialised equipment and expertise, and are not accessible to every biologist. PALM/STORM, STED and 3D-SIM require image reconstruction/processing as a final step and can be applied in a restricted manner to live imaging [START_REF] Schermelleh | A guide to super-resolution fluorescence microscopy[END_REF]. However, even though these techniques have been democratised recently by the emergence of commercial solutions they ask for high technical expertise and specific hardware, they are time consuming and display inherently high photo-toxicity [START_REF] Schermelleh | A guide to super-resolution fluorescence microscopy[END_REF].

On the other hand, confocal microscopy, a valuable tool for optical sectioning of thick samples has become a standard technique accessible to every biologist. A conventional confocal microscope may reach resolutions of 200 nm laterally and 500 nm axially, if properly adjusted [START_REF] Pawley | Fundamental Limits of Confocal Microscopy, Handbook of Biological Confocal Microscopy[END_REF]. Furthermore, image restoration by deconvolution, a powerful computational process used to reduce out-of-focus light in 3D-fluorescent images, can be combined with confocal imaging. This combination has been proven useful to enhance image contrast, improve the signal-to-noise ratio and resolution [START_REF] Cannell | Image enhancement by deconvolution[END_REF][START_REF] Shaw | Comparison of widefield/deconvolution and Confocal microscopy for Three-Dimensional Imaging[END_REF][START_REF] Heck | A deconvolution method to improve automated 3D-analysis of dendritic spines: application to a mouse model of Huntington's disease[END_REF].

However, deconvolution of confocal data is still scarcely used. This might be due to the fact that deconvolution requires expertise and a neat measurement of the impulse response of the optical system; the point spread function (PSF). When focussing deep into a biological specimen the PSF is distorted by spherical aberrations due to refractive index mismatches [START_REF] Egner | Aberrations in confocal and Multi-photon Fluorescence Microscopy induced by refractive index mismatch[END_REF], and might then not be precise enough to perform high fidelity deconvolution. Besides, scattering is also an issue in thick and optically dense samples. To counteract the negative influence of the distorted PSF on the deconvolution results, one might use a PSF taken directly from the biological image, if the structures investigated have the appropriate shape [START_REF] Von Tiedemann | Image adaptive point-spread function estimation and deconvolution for in vivo confocal microscopy[END_REF]. It is also possible to compute theoretical PSF [START_REF] Dieterlen | Validation of image processing tools for 3-D fluorescence microscopy[END_REF].

Alternatively, adaptive optics could be used to minimize the distortion of the wave front in depth of the sample [START_REF] Booth | Aberrations and adaptive optics in super-resolution microscopy[END_REF]. Another possibility to avoid spherical aberration and thus distortion of the PSF would be refractive index matching of the sample [START_REF] Staudt | 2,2'-thiodiethanol: a new water-soluble mounting medium for high-resolution optical microscopy[END_REF]. We have recently shown that accurate refractive index matching of the sample and the mounting medium improves considerably the axial resolution and depth penetration in fixed brain tissue [17]. We furthermore showed that the fixed biological tissue has an overall refractive index that is close to that of the optical system (1.518), rendering the tissue very transparent. This minimizes scattering and the PSF measured after focussing up to 120 µm into this sample was not impaired by spherical aberration and thus distortion-free.

In the present study, we aim at applying this knowledge on the mouse oocyte, a rather bulky sample with a diameter of nearly 100 μm. Our scientific project aims at improving microtubule (MT) imaging at high resolution in the meiotic spindle. High resolution imaging of the meiotic spindle will be important to tackle outstanding questions regarding the influence of the structural organisation of the MTs network on meiotic division in the future. Microtubules of the meiotic spindle are hollow tubes of 25 nm in diameter, consisting of 13 protofilaments composed of alpha-and beta-tubulin heterodimers. These polymers are highly dynamic and stochastically switch between polymerisation and depolymerisation, a process known as dynamic instability [START_REF] Kirschner | Microtubule dynamics[END_REF]. MT dynamics and structural organisation is tightly regulated by a plethora of MT-associated proteins (MAPs) [START_REF] Alfaro-Aco | Building the Microtubule Cytoskeleton Piece by Piece[END_REF]. To understand how the precise organisation of the MT network contributes to essential cellular processes such as chromosome segregation during meiosis, a detailed reconstruction and analysis of the spatial distribution of MTs as well as their length, orientation and anchoring to kinetochores in oocytes is decisive. However, conventional fluorescence microscopy of these polymers and their associated proteins in mammalian oocytes is limited by light diffraction and by the thickness of the sample.

Visualizing the attachment of individual microtubule fibres on each sister kinetochore is beyond the resolution limit of classical confocal imaging.

In the present paper, we combine refractive index matching, optimisation of imaging parameters and deconvolution of confocal data. Several commercial solutions have been recently developed to improve lateral and axial resolution of confocal imaging. Leica's Hyvolution and Olympus' FV-OSR solution both use a standard confocal microscope, customized with optimized acquisition and detection system (pinhole closure and high sensitivity, low-noise detector) coupled to image restoration. A different approach is used in the Zeiss' Airyscan configuration where all photons are collected from the volume excited with concentrically arranged hexagonal detector array. Its detection area consists of 32 single detector elements, acting as very small pinholes. The signals from these detector elements are reassigned to their correct position, producing an image with increased signal-to-noise ratio and resolution.

We propose a workflow similar to the solution from Leica and Olympus, with the difference, that we adapt the mounting medium by refractive index matching with the biological sample. This matching decreases the impact of the spherical aberration and improves thus axial resolution. Our workflow allows to improve the lateral and axial resolution of the confocal microscope, gaining resolutions comparable to 3D-SIM-imaging, but with a much better depth penetration. In a first step, we characterize the resolution of our optical system by acquiring sub resolution fluorescent beads using this workflow. We then demonstrate that this optimized workflow improves considerably the resolution of biological imaging when applied to the meiotic spindle of mouse oocytes.

Materials and Methods

Oocyte culture and whole-mount immunofluorescence

Oocytes were obtained from ovaries of adult Swiss mice 10 to 16 weeks old (Janvier, France) and cultured in self-made M2 medium, as previously described [START_REF] Touati | Mouse oocytes depend on BubR1 for proper chromosome segregation but not for prophase I arrest[END_REF]. Oocytes entering meiosis I, visible through Germinal Vesicle Breakdown (GVBD), within 60 min after harvesting were used for experiments. Before using oocytes for fixation, the zona pellucida was removed by incubation in Tyrode's solution [START_REF] Chambon | Chromosome spreads with centromere staining in mouse oocytes[END_REF].

Between 6h and 6h30 after GVBD, zona pellucida-free oocytes were placed into self-made chambers in order to manipulate them easily. Oocytes were cytospinned in chambers coated with concanavalin A (at 0.2mg/ml, from Sigma) to keep oocytes in place, for 13 min at 1400 RPM at 38°C. Cold-stable spindles were obtained by incubating the individual chambers on an ice-water bath for 4 min before fixation. Oocytes were fixed with 2 % Formaldehyde solution (Sigma-Aldrich; F1635) in BRB80 buffer with Triton TX-100 during 30min at 38°C, as described [START_REF] Hached | Mps1 at kinetochores is essential for female mouse meiosis I[END_REF]. After PBS washing, oocytes were incubated overnight in a PBS-BSA 3% Triton TX-100 solution for permeabilisation and blocking unspecific antibody binding.

Antibody staining of fixed oocytes was done using the antibodies below at the indicated dilutions in PBS: Human CREST autoimmune antibody (Cellon SA, HCT-0100; 1/100), and mouse monoclonal αtubulin (DM1A) coupled to FITC (Sigma-Aldrich, F2168; 1/100). As secondary antibody CY3 anti rabbit (Jackson ImmunoResearch, 711-166-152; 1/200) was used. Antibody incubation times were at least one hour. Chromosomes were stained with 5 µg/ml Hoechst 33342 during incubation with the secondary antibody. Oocytes were washed several times with PBS between antibody incubations and before mounting the chamber [START_REF] Hached | Mps1 at kinetochores is essential for female mouse meiosis I[END_REF].

Mounting media

We used AF1 (Citifluor, UK), a commercially available mounting medium with a refractive index of 1.463 and AF1+, a modified AF1 solution harbouring a refractive index of 1.518. The refractive index increase of AF1+ solution was obtained by adding 83 % (w/w) of Methyl-Phenylsulfoxid (Sigma-Aldrich, #261696) to AF1-solution. Refractive indices were verified at 21°C using a refractometer (Mettler Toledo, Switzerland).

Preparation of Sub-resolution fluorescent beads

The experimental setup for evaluating lateral and axial resolution at the coverslip and in depth was designed as followed: fluorescent beads (PS-Speck, Lifetechnologies), of a diameter of 100 nm, and loaded with yellow-green fluorescent dye, were diluted in water (1/800 v/v). Drops of the water-diluted sample were put on the surface of the coverslip (Menzel Glaeser #1.5, Agar scientific) or slide and air-dried. The coverslip and slide were separated by one layer of adhesive tape (Scotch®, 3M) with a nominal thickness of 60 µm and the volume filled with a drop of the respective mounting medium.

Confocal laser scanning microscopy

8-bit Images were collected using a Leica 63x oil immersion objective (HCX Plan APO CS, NA 1.4, working distance 0.14 mm) with an inverted Leica laser-scanning confocal microscope TCS SP5 II (Leica Microsystems, Heidelberg, Germany) equipped with a GaAsP hybrid detection system at a sampling rate of 40 nm in x,y-and 83 nm in z-direction unless otherwise stated. Fluorochromes were detected using laser lines 405 nm, 488 nm and 561 nm. Imaging was performed in a temperaturecontrolled room at 21°C. The detection pinhole aperture was adjusted to 1 Airy unit or 0.6 Airy units.

Slow scanning speeds may improve the signal-to-noise ratio and the resolution of the image since more photons may be collected by pixel. The cutoff value has to be determined for every biological sample, since slower scan speed results in higher photo bleaching. We tested scan speeds of 400, 200, 100 and 50 Hz. For the given biological sample the optimal scan speed was 400 Hz. We also tested Immersion oils with refractive indexes from 1.510 to 1.520 to minimize spherical aberration (supp. figure 2). We used an immersion oil harbouring a refractive index of 1.518 for all image acquisitions.

Image acquisition of beads

Bead images were obtained as in [17] with the following modifications. GaAsP gain was set to 15% and the laser power adjusted so that the signal occupied the full dynamic range of the detector, but saturated voxels were carefully avoided. Beads were imaged starting and finishing the stack at least 5μm below and above the bead centre. Beads were visually checked and improper stacks were discarded before determining the microscope PSF.

Beads used for resolution measurements (Fig. 1) were acquired at a sampling rate of 32 nm in x, y and 42 nm in z-direction. An average of 60 beads was taken for each measurement. Beads used for deconvolution were acquired at a sampling rate of 40 nm in x,y and 83 nm in z-direction. At least 15 beads were registered and averaged in order to increase the SNR for deconvolution.

Resolution measurement with beads

Resolution measurements were carried out as in [START_REF] Kirschner | Microtubule dynamics[END_REF] using ImageJ with the following modifications:

We took the maximum intensity plane of the image stack along the x, y and z-axis resulting in 1D intensity profiles. A x,y-and z-line-profile was generated manually along a line passing through the 2D maximum. From this image, fluorescence intensity profiles were fitted to a Gaussian curve, using ImageJ's built-in curve fitting function. The full width at half maximum (FWHM) of the gaussian curve was calculated for each profile, based on the parameters retrieved from the fitting (supp. fig. 3 A).

Deconvolution and image treatment:

Confocal images of beads and biological data were deconvoluted with the Huygens 3.7 software (Scientific Volume Imaging, Hilversum, Netherlands) using a measured PSF (supp. fig. 3B). We tested the Classical Maximum Likelihood Estimation (CMLE) algorithm and the Good's (GMLE) algorithm.

Signal-to-noise ratios, also called the R-parameter from 15 to 20 are recommended by the manufacturer for noisy confocal images R-values > 20 for low noise wide-field images. Since we used a confocal microscope with a low noise detector, we tested R-parameters of 26 to 34 and visually inspected the results. In images of figure 2 to 6 brightness and contrast were adjusted equally for all images after deconvolution and before 3D reconstruction and volume rendering. Image analysis was performed with ImageJ [23] and 3D-volume rendering was performed with ICY [START_REF] De Chaumont | Icy: an open bioimage informatics platform for extended reproducible research[END_REF]. Deconvolution was perfomed with a measured PSF, with a quality threshold of 0.001, varying R-parameters as stated, performing 500 iterations in the optimized mode with CMLE-method and 125 iterations with the GMLE-method. Denoising was perfomed using the Pure Denoise Plugin [START_REF] Luisier | Fast Interscale Wavelet Denoising of Poisson-corrupted Images[END_REF]26,[START_REF] Luisier | Fast Haar-Wavelet Denoising of Multidimensional Fluorescence Microscopy Data[END_REF] implemented in ImageJ with the cycle spin and the multiframe parameter set to maximum.

Signal to Noise Ratio calculation

To calculate the SNR, we choose a threshold with the Shanbhag mode on ImageJ (supp. figure 4). We created a mask to measure the average intensity of our structure and with the inverted thresholded image; we create a second mask to measure the average intensities of the background. We finally calculate the SNR as average intensity/average background.

Results and discussion

Resolution measurement using sub-resolution microspheres

We first aimed at developing a workflow, consisting of optimized sample preparation, image acquisition and deconvolution parameters to obtain the best resolution with our standard confocal microscope. We therefore investigated the lateral and axial resolution of the confocal microscope, using a high numerical aperture lens (63x, NA 1.4, oil immersion).

We determined the Point Spread Function (PSF) using two different mounting media, and varying the detection pinhole aperture. The PSF is the impulse response of the focused optical system and gives information about the lateral and axial resolution of the optical system. The PSF was measured by imaging 3D stacks of 100 nm fluorescent beads at different depths. The lateral and axial resolutions were estimated from the full-width at half-maximum (FWHM) of the intensity profile of the PSF. We compared beads mounted in AF1 with a refractive index of 1.463 with beads mounted in AF1+, a modified AF1, harbouring a refractive index of 1.518. In the first case we introduced a mismatch of the refractive index (ri) between the sample (ri=1.463) and the glass-oil-interface (ri=1.518) thus producing spherical aberrations when focussing away from the coverslip. In the second case, we perfectly matched the sample (ri=1.518) with the glass-oil-interface (ri=1.518), thus minimizing spherical aberrations.

Resolution improvement by deconvolution

Firstly, we determined the resolution of the confocal microscope before and after deconvolution of confocal data of beads in the ri-mismatched and ri-matched case with the detector pinhole set at 1 Airy Unit (AU). The resolution of the optical system taken from the raw images was 225 nm ± 12 nm in x,y and 546 nm ± 18 nm in z-direction in AF1-medium and 223 nm ± 11 nm in x,y and 550 nm ± 19 nm in z-direction in AF1+-medium (Fig. 1A,table 1). After deconvolution, we observed a 1.7-fold increase in lateral resolution for both mounting media, a 2.8-fold increase in axial resolution for AF1 and a 2.5-fold increase for AF1+ for beads imaged next to the coverslip (Fig. 1 A, table 1). The resolution values for AF1 and AF1+-mounted beads were comparable, corroborating previous findings, that close to the coverslip, a refractive index mismatch does not impair lateral and axial resolution since high numerical aperture objectives are corrected for this [START_REF] Von Tiedemann | Image adaptive point-spread function estimation and deconvolution for in vivo confocal microscopy[END_REF][START_REF] Kirschner | Microtubule dynamics[END_REF]. Having established that deconvolution improves indeed the resolution of confocal data, we secondly analysed how the refractive index mismatch would impact resolution when looking at beads at a depth of 60 μm (Fig 1B, deconvoluted data, table 1). The lateral resolution of the beads was comparable for the two mounting media, which confirmed that the lateral resolution in depth is not impaired in a mismatched system [START_REF] Egner | Aberrations in confocal and Multi-photon Fluorescence Microscopy induced by refractive index mismatch[END_REF]17]. However, axial resolution was 2.2-times inferior in AF1-medium compared to AF1+-medium. Besides, we observed a considerable loss of peak intensity at a depth of 60 μm in AF1-medium. We had to adjust the laser power accordingly in order to exploit the full dynamic range of the image. This confirms the data, we published recently for beads measured in CFM3, another high refractive index mounting medium (ri=1.518) where we observed the same phenomenon [17]. The axial resolution of beads remained constant in depth compared to the coverslip when applying refractive index matching with AF1+. This means that spherical aberration is minimized when using AF1+ medium, and our observation corroborated previous theoretical and experimental findings of Hell and co-workers [START_REF] Egner | Aberrations in confocal and Multi-photon Fluorescence Microscopy induced by refractive index mismatch[END_REF] and our own findings for a perfectly matched system [17].

Resolution improvement by closing the pinhole aperture

We then wanted to test if we could increase resolution of the confocal microscope by closing the detection pinhole. We compared resolution at the coverslip and in a depth of 60 μm with the detector pinhole set to 0.6 AU (Fig. 1C,table 1) and after deconvolution of the data. Our choice of the 0.6 AU pinhole size was based on several tests on our biological data. We acquired the same type of biological sample with different pinhole sizes, from 1 AU to 0,4 AU and observed that 0.6 AU is the threshold where we discard enough diffraction signal without photo-bleaching and with a good contrast. Since the result depends largely on the quality and the photo-stability of the biological sample the optimal pinhole value has to be evaluated for each biological sample.

We observed lateral resolutions of 120 nm ± 5 nm (coverslip) and 136 nm ± 9 nm (depth) and axial resolutions of 182 nm ± 10 nm (coverslip) and 471 nm ± 89 nm (depth) for the AF1 medium (Fig. 1C, table 1). We had again to compensate for the loss of peak intensity by adjusting the laser power with the AF1 medium as described before. For the AF1+-medium, lateral resolutions of 106 nm ± 4 nm (coverslip) and 120 nm± 7 nm (depth) and axial resolutions of 164 nm ± 8 nm (coverslip) and 192 nm ± 17 nm (depth) were measured. Fig. 1D andE show the direct comparison of the optical resolution at 1 AU and 0.6 AU close to the coverslip (Fig. 1D, table 1) and in 60 μm depth (Fig. 1E, table 1), respectively. We observed a 1.3-fold increase in lateral and axial resolution measured in AF1-medium close to the coverslip (Fig. 1D, table 1) when setting the detection pinhole to 0.6 AU. The use of AF1+ medium improved resolution 1.3fold laterally and 1.4-fold axially at AU 0.6 at the coverslip (Fig1. D, table 1). In depth, a significant lateral and axial resolution improvement was observed only for AF1+-medium at 0.6 AU. Besides, the use of AF1+ allowed maintaining a good axial resolution of 192 nm ± 17 nm that is 2.5-times superior than with AF1-medium (471 nm ± 89 nm), when closing the pinhole. It is thus advantageous to use refractive index matching when closing the detector aperture and imaging in depth.

Closing the pinhole indeed increased lateral and axial resolution 1.3 to 1.4-fold. These results are in good agreement with data measured by Cox and Sheppard [START_REF] Cox | Practical limits of resolution in confocal and non-linear microscopy[END_REF], who observed a 1.4-fold increase in resolution after closing the pinhole aperture to 0.5 using a Leica TSC SP2 confocal microscope.

By optimizing sample preparation, image acquisition parameters and performing deconvolution, our workflow allowed us to obtain a considerable gain in lateral and axial resolution throughout the sample thickness. Firstly, the lateral resolution attained with 100 nm sub resolution fluorescent beads close to the coverslip attained a value of 106 nm ± 4 nm. This is close to the real size of the beads and comparable to the values published earlier for 3D-SIM-imaging [START_REF] Gustafsson | Three-dimensional resolution doubling in wide-field fluorescence microscopy by structured illumination[END_REF]. Secondly, the axial resolution we obtained with our workflow was 164 nm ± 8 nm and is thus approximately 1.6-fold better than the values published for 3D-SIM by Gustafsson [START_REF] Gustafsson | Three-dimensional resolution doubling in wide-field fluorescence microscopy by structured illumination[END_REF]. However, these results are in good agreement with a paper by Schrader and colleagues [START_REF] Schrader | Potential of confocal microscopes to resolve in the 50-100 nm range[END_REF], who demonstrated earlier that a resolution of 80 nm could be obtained axially on 50 nm gold particles with a confocal microscope and after deconvolution. Given that we used 100 nm beads, our results are consistent with these previous findings. Thirdly, refractive index matching of the sample allowed to maintain 3D-SIM-like resolution of 120 nm ± 7 nm laterally and 192 ± 17 nm axially in a depth of 60 μm. This depth surpasses by a factor 3 the imaging depth of standard 3D-SIM [6].

Deconvolution optimization in biological data

Having established optimal image acquisition parameters, we then wanted to apply our workflow on biological samples and chose the meiotic spindle of the mouse oocyte labelled with alpha-tubulin coupled with FITC. We acquired 3D-images of meiotic spindles with a pinhole size of 0.6 Airy units at a sampling rate of 40 nm in x,y and 83 nm in z-direction. We wanted to determine the optimal deconvolution parameters and tested the Classic Maximum Likelihood Estimation (CMLE) and the Good's roughness Maximum Likelyhood Estimation algorithm (GMLE), implemented in Huygens deconvolution software. CMLE and GMLE are iterative restoration method optimizing the likelihood of an estimate of the object given the measured image and the PSF. We also studied the influence of denoising prior to deconvolution using the Pure Denoising plugin implemented in ImageJ [26,[START_REF] Luisier | Fast Haar-Wavelet Denoising of Multidimensional Fluorescence Microscopy Data[END_REF][START_REF] Cox | Practical limits of resolution in confocal and non-linear microscopy[END_REF].

Finally we compared raw data and data optimized by our workflow after 3D-reconstruction.

Establishment of the optimal Regularisation Parameter (R-parameter)

The Huygens software uses what the manufacturer calls the Signal-to-Noise ratio as a Regularization Parameter (R-parameter), i.e. as a parameter that controls the sharpness of the restoration result.

We will employ the term R-parameter in the following in order to avoid confusion with the image intrinsic Signal-to-Noise ratio (SNR) that we calculated to evaluate image quality.

The R-parameter is calculated as the square root of the number of photons in the brightest part of the image and controls the sharpness of the restoration result. The higher the R-parameter, the sharper is the deconvolution result. R-parameter values recommended by the manufacturer for noise-prone confocal data lie between 15 and 20. Since we used a low-level-noise detector for confocal imaging, we tested R-parameter values between 16 and 34. After adjusting the Rparameter, we visually inspected the deconvolution result and observed that at an R-parameter >26, microtubules were separated optimally, when compared to the raw image (Fig. 2 A-D). We also carried out fluorescence intensity profile measurements to obtain a measurement of the separation of tubulin structures in a single plane (Fig. 2A'-D4) and calculated the image inherent SNR before and after deconvolution to estimate image improvement. We noticed that microtubule separation attained a maximum, setting the R-parameter to 30 (Fig. 2 A'-D'). The SNR increased from 3,1 in the raw image to 3,32 after deconvolution with R-parameter set to 26. We started separating 5 peaks in the region of interest at an R-parameter of 28 (Fig. 2C',box). At higher R-parameter, the SNR decreased slightly to 3.24 (R-parameter 28) and 3.20 (R-parameter [START_REF] Boutet De Monvel | Image restoration for confocal microscopy: improving the limits of deconvolution, with application to the visualization of the mammalian hearing organ[END_REF], showing that the Rparameter introduces sharpness into the image. When using R-parameters >30 we started introducing deconvolution artefacts such as patterning, indicating that an R-parameter of 30 is optimal for our dataset (Fig 1, Supplementary Data 1 A-C).

Influence of denoising on the deconvolution of confocal data

Having established the optimal R-parameter, we now aimed at studying the impact of denoising before deconvolution. It has been proposed that denoising may be crucial before deconvolution and may improve the restoration result [START_REF] Boutet De Monvel | Image restoration for confocal microscopy: improving the limits of deconvolution, with application to the visualization of the mammalian hearing organ[END_REF]. We therefore denoised the confocal data with the "Pure Denoising" plugin implemented in ImageJ [START_REF] Luisier | Fast Interscale Wavelet Denoising of Poisson-corrupted Images[END_REF]26,[START_REF] Luisier | Fast Haar-Wavelet Denoising of Multidimensional Fluorescence Microscopy Data[END_REF]. This plugin is well adapted for Poisson noise, which is characteristic of fluorescence images. It estimates automatically noise parameters such as detector gain, detector offset and the standard deviation of the additive white Gaussian noise. We adjusted the number of cycle-spins and the number of adjacent frames used to maximum values in order to get the highest quality for our images. Denoising with the quality criteria set to maximum was rather computing intense and took several hours for a 3D-image stack. Fig. 3 shows a comparison of the raw image (A, A'), the denoised image (B, B'), the deconvolution of the raw image (C, C') and the deconvolution of the denoised image (D, D'). When inspecting the fluorescence intensity profiles, denoising alone did not improve the resolution of the image (Fig. 3B,B'), even if it improved considerably the SNR. We observed a smoothing of the fluorescence intensity profile (Fig. 3B',box). We then compared the deconvoluted dataset (Fig. 3C,C',box) with the denoised and deconvoluted image (Fig. 3D,D',box). We observed 5 main peaks in the region of interest of the deconvoluted image (Fig. 3C'). Denoising improved SNR of the restoration result. After denoising we were still able to distinguish 5 peaks (Fig. 3D, D' and box), even if the indentation was less prominent than in the non-denoised dataset (Fig. 3C,C',box). The overall aspect of the image was smoother and as a result microtubules appeared less spotty after denoising (Fig. 3 D, insert).

Comparison of CMLE and GMLE algorithm for the deconvolution of confocal data

In a last step, we wanted to compare the CMLE-algorithm with the Goods Maximum Likelihood Estimation (GMLE) algorithm. The GMLE-methods needs 4-times less iterations than the CMLEmethod to obtain comparable restoration results, thus saving computation time. We compared the image after deconvolution with the CMLE-algorithm (Fig. 4A,A') with that deconvoluted with the GMLE-algorithm (Fig. 4B,B'). After visual inspection, we observed that the GMLE-algorithm smoothed the deconvolution result more that the CMLE-algorithm. Both algorithms gave comparable results in separating the microtubules and 5 main peaks could be distinguished in both cases (Fig. 4A', B', box), however, the SNR was higher for the GMLE-algorithm. Denoising did improve the SNR of the image deconvoluted with CMLE-algorithm, however, the SNR after GMLE-deconvolution remained constant, indicating that denoising did not further improve image quality in this situation (Fig. 4,D,D').

Influence of refractive index matching on the 3D-image resolution

We have shown previously, that refractive index matching greatly improves axial resolution of confocal data in the mouse brain [17]. We therefore mounted mouse oocytes in AF1 and AF1+, respectively and performed 3D imaging. Fig. 5 shows lateral and axial views of the tubulin staining and after deconvolution. We observed a 2-fold increase of lateral resolution when comparing raw data to deconvoluted data. When considering the lateral maximum projection, the overall aspect of the image in the two mounting media was comparable (Fig. 5 A,B). We then measured the diameter of microtubule bundles throughout the depth of the image stacks and detected comparable lateral resolution for AF1 and AF1+ (Fig. 5 A,A'). However, axial resolution seemed to be better in the rimatched sample with AF1+ after visual inspection (Fig. 5 B'). We measured axial resolutions and found indeed a 1.4-fold increase of resolution in AF1+-medium, meaning that refractive index matching is crucial in the biological sample. This confirms the data obtained with beads, where a comparable gain of axial resolution has been observed in depth using the ri-matched system (Fig. 1C) and our data on biological samples published earlier [17].

3D-reconstruction of the meiotic spindle

Having established all optimal parameters, we applied them to acquire 3D-stacks of triple stained mouse oocytes. We evaluated the overall aspect of the meiotic spindle resolution. We used alphatubulin coupled to FITC, to mark microtubules, Crest, a protein staining the kinetochore, which was revealed with a secondary antibody coupled to Cy3, and Hoechst to stain the chromosomes. Fig. 6 shows the 3D rendering of the raw data (A) and the deconvoluted data (B) of the meiotic spindle. We observed a clear improvement in resolution and much better separated microtubules (Fig. 6 B', arrows) and the kinetochores (Fig. 6 B', arrowheads) as compared to the raw data (Fig. 6 A',, arrows and arrowheads). With these results, we provide evidence that our workflow, consisting of refractive index matching, optimized image acquisition and deconvolution improves greatly the overall resolution of the meiotic spindle image.

Conclusion

The proposed workflow provides an important benefit to confocal imaging of fixed samples with a high numerical aperture lens. We show first with beads that deconvolution improves by 1.7-fold lateral and by >2.6-fold the axial resolution of confocal data. Secondly, we provide evidence that refractive index matching is particular important to improve resolution and signal strength of the confocal system, when imaging in depth, where we obtain a 2.2-fold axial resolution improvement for the matched system. Thirdly, we show that closing the pinhole improves lateral and axial resolution by a factor 1.3 close to the coverslip and that the spatial invariance of the PSF in depth is crucial to preserve resolution in depth. We observe that deconvolution of the biological confocal data improves 2-fold the resolution compared to raw confocal data and that refractive index matching helps to maintain optimal axial resolution in the biological sample. We lastly show that denoising is advantageous for image quality when using the CMLE-algorithm, however, it is time consuming. Denoising is not necessary when using the GLME-algorithm, which also needs much less computing time. Summarized from these results, we were able to obtain confocal images having a 3D-SIM-like resolution in a rather bulky sample, the mouse oocyte. It would be of further interest to apply our workflow on other biological samples and see if we are able to maintain a 3D-SIM-like resolution at the full working distance of the objective (140 μm).

Acknowledgements

The authors thank Jean-François Gilles and Nicolas Heck for helpful discussions and Marc Dos Santos for precious help with statistical analysis and Adobe Illustrator. Work on mouse oocytes was supported by an ANR Grant (ANR-12-BSV2-0005-01) to KW. Confocal work was carried out at the Institute of Biology Paris-Seine Imaging Facility that is significantly (or: which receives significant support from...) supported by the "Conseil Regional Ile-de France", the French national research council (CNRS) and Sorbonne University, UPMC Univ Paris 06.

Figure legends:

Fig. 1: Resolution measurements with sub-resolution microspheres: influence of deconvolution, refractive index and pinhole aperture 100 nm green fluorescent beads were used to measure the lateral and axial resolution of the optical system at the coverslip and in 60 µm depth. For excitation, the 488 nm laser line was used and two series with the detection pinhole set to 1 Airy Unit (AU) and 0.6 AU were acquired. The Full Width at Half Maximum (FWHM) was measured using a line profile in x, y-direction and x, z-direction in the maximum intensity plane of the image stack, respectively and transferred to a graph (A-E). All data in brackets are given as mean ± SD and have the unit nm.

(

  A) Raw bead images (black and dotted bars) and deconvoluted bead images (striped and white bars) taken at the coverslip were compared. Deconvolution improved significantly the lateral (raw: 225 ± 12; decon: 134 ± 3) and axial resolution (raw: 546 ± 18; decon: 193 ± 13) of beads in AF1. Deconvolution improved also significantly the lateral (raw: 223 ± 11; decon: 134 ± 3) and axial resolution (raw: 550 ± 19; decon: 224 ± 18) of beads in AF1+-medium. (B, C) Lateral and axial resolution of deconvoluted bead images at the coverslip (black and dotted bars) and in 60 µm depth (striped and white bars) at 1 AU (B) and 0.6 AU (C). Resolution measurements in depth show a significant axial resolution loss in AF1-medium at 1 AU (coverslip: 193 ± 13; depth: 453 ± 89) and at 0.6 AU (coverslip: 182 ± 10; depth: 471 ± 89). (D, E) Comparison of the resolution of beads taken at 1 Airy unit (black and dotted bars) and 0.6 Airy unit (striped and white bars) at the coverslip (D) and in 60 µm depth (E). Resolution measurements at the coverslip (D) show a significant lateral resolution improvement when closing the pinhole to 0,6 AU in AF1-medium (1 AU: 134 ± 3; 0.6 AU: 120 ± 5, Mann-Whitey test, p= 0.0001), and a substantial axial resolution increase (1 AU: 193 ± 13; 0.6 AU: 182 ± 10, Mann-Whitey test, p= 0.0001). A significant lateral resolution improvement at 0.6 AU is also observed in AF1+-medium (1 AU: 134 ± 3; 0.6 AU: 106 ± 4, Mann-Whitey test, p= 0.0001) and axial resolution improves substantially, too (1 AU: 210 ± 21; 0.6 AU: 192 ± 17, Mann-Whitey test, p= 0.0001).In depth (E), no resolution improvement is observed laterally (1 AU: 140 ± 10; 0.6 AU: 136 ± 9) or axially (1 AU: 453 ± 89; 0.6 AU: 471 ± 89) in AF1-medium. However, when using AF1+-medium, resolution improves significantly in the lateral direction (1 AU: 133 ± 6; 0.6 AU: 120 ± 7, Mann-Whitey test, p= 0.0001) and in the axial direction (1 AU: 210 ± 21; 0.6 AU: 192 ± 17, Mann-Whitey test, p= 0.0003).

Fig. 2 :

 2 Fig. 2: Determination of optimal deconvolution parameters, R-Parameter Mouse oocytes stained with FITC-tubulin and imaged at 488 nm with the detection pinhole adjusted to 0.6 Airy Unit and taken with a 40 nm pixel size and an 83 nm z step size were deconvoluted varying the R-parameter between 26, 28 and 30. Images show x, y-maximum projections of the raw data (A), and the deconvoluted data with R-parameter 26 (B), R-parameter 28 (C), and R-parameter 30 (D). Inserts show the position of the line profile used for fluorescence intensity measurements shown in the graphs next to each image (A'-D'). S/N-ratios indicated in (A'-D') are a value for the image intrinsic noise and were calculated as described in (2.8). The boxes indicate the regions of interest, where 5 main peaks can be separated after deconvolution. Note that R-parameter of 30 gave best results. Scale bar is 5 µm.

Fig. 3 :

 3 Fig. 3: Influence of denoising on the deconvolution result Images show x, y-maximum projections of the raw data (A), the deconvoluted data (B), the denoised data (C) and the denoised, deconvoluted data (D). Deconvolution was carried using the CMLEalgorithm (R-parameter=30, 500 iterations). Inserts show the position of the line profile used for

Fig. 4 :

 4 Fig. 4: Comparison of deconvolution algorithm CMLE and GMLE Images show data deconvoluted with two different deconvolution algorithms implemented in the Huygens software, the Classical Maximum Likelyhood Estimation (CMLE) and the Good's Maximum Likelyhood Estimation (GMLE). Deconvolution was carried with the R-parameter set to 30 for both methods. Images show x, y-maximum projections of the CMLE-deconvoluted data (A), the GMLE deconvoluted data (B) the denoised, CMLE deconvoluted data (C), and the denoised, GMLE deconvoluted data (D). Inserts show the position of the line profile used for fluorescence intensity measurements shown in the graphs next to each image (A'-D'), the boxes indicate the position of the 5 main peaks compaired before. Note that denoising improves the S/N of the CMLE-method, but not of the GMLE-method. The GMLE-algorithm gives the best visual result. Scale bar is 5 µm.

Fig. 5 :

 5 Fig. 5: Impact of refractive index matching in mouse oocytesImages show x, y-maximum projections (A, B, lateral view), a cut-out of a single plane in a depth of 10 µm (insert of A, B) and x, z-projections (A', B') of mouse oocytes stained with FITC-tubulin and imaged at 488 nm in AF1-medium (A, A') and AF1+-medium (B, B'). All data in brackets are given as mean ± SD and have the unit nm.Microtubule bundle diameters were measured in single images of raw data (images not shown) and deconvoluted data (raw data/ deconvolved data; 390 ± 131 /218 ± 45). A 2-fold amelioration of resolution was noted after deconvolution. Lateral resolution of microtubule bundles in AF1 and AF1+-medium was comparable (AF1/ AF1+; 232 ± 33/204 ± 57). Axial resolution of microtubule bundles was measured in AF1 and AF1+-medium (AF1/ AF1+; 402 ± 107 /292 ± 49). Note that refractive index matching with AF1+ improves considerably the axial resolution (A', B'). Scale bars are 5 µm.

Fig. 6 :

 6 Fig. 6: 3D-reconstruction of data Mouse oocytes were triple stained with Hoechst (chromosomes), FITC (microtubules) and CY3 (kinetochores). Images show 3D-renderings of raw data (A-A') and deconvoluted data (B-B'). Magnifications of a region of (A, B) with arrows indicating microtubules (A', B') and arrowheads

Supplementary data 1:

Images show data deconvoluted with two different deconvolution algorithms implemented in the Huygens software, the Classical Maximum Likelyhood Estimation (CMLE) and the Good's Maximum Likelyhood Estimation (GMLE). Deconvolution was carried with the R-parameter set to 30, 32 and 34 for both methods. Images show x, y-maximum projections of the CMLE-deconvoluted data (A) and the GMLE deconvoluted data (B). Inserts show magnifications of microtubules. R-parameters >30 improve the S/N of the GMLE-method, but not of the CMLE-method. When using R-parameters >30 we started introducing deconvolution artefacts such as patterning when using the CMLE-method (inserts), indicating that an R-parameter of 30 is optimal for our dataset. Scale bar is 5 µm. Supplementary data 2: Influence of immersion oil refractive index on axial resolution of subresolution fluorescent beads. Note that axial resolution is best with immersion oil harbouring refractive index of 1.516-1.518.

Supplementary data 3: Determination of bead resolution (A) and xz view of beads (B) comparing

their PSF when images at the coverslip and at depth in both non-matched (AF1) and matched RI (AF1+).

Supplementary data 4: To calculate the SNR, images were thresholded with the Shanbhag mode on ImageJ. We create a mask from the duplicated image. Two selections are necessary: one for measuring the average intensity of the signal, the other for measuring the background. We finally calculate the SNR as average intensity signal/average intensity background.