Batoul Dahman
email: batoul.dahman@univ-rennes1.fr

Francis Bessière

Jean-Louis Dillenseger

Ultrasound to CT rigid image registration using CNN for the HIFU treatment of heart arrhythmias

Keywords: Rigid-body image registration, Ultrasound Imaging-cardiac, Multimodal image fusion

Image-guided thermal ablations have become an important therapeutic option for patient with cardiac arrhythmia, it is minimally-invasive and provides better and faster patient recovery. However, to enhance the ablation guidance, the therapist needs to link by image registration the intraoperative images to the high-resolution anatomical preoperative imaging, in which the ablation path has been defined. In this work, we present a convolutional neural networks (CNNs) framework for transesophageal ultrasound/computed tomography image registration to solve the problem of high computation time of the classical iterative methods, which is not suitable for a real-time application. We propose the following process: we first pass the input moving and fixed image pairs through a siamese architecture consisting of convolutional layers, thus extracting features of moving and fixed maps analogous to dense local descriptors, then matching the feature maps, and finally pass this correspondence feature map into a registration network, which directly outputs the registration parameters set of the rigid registration. Accuracy of the registration is quantified based on the Target Registration Error (TRE) for specific anatomical landmarks. Results of the registration process show a median TRE of 2.2 mm for all the fiducial points, and the registration computation time was around 3 ms comparing to the classic iterative methods which takes around 70 seconds for one image pair. In our future work we are going to perform our approach on 2D/3D learning-based registration to refine the estimation of the transesophageal probe pose in the 3D preoperative volume.

INTRODUCTION

In the last 20 years, cardiac arrhythmias (heart rhythm problems) have become one of the most important public health problems and a significant cause of increasing health care costs in Western countries. Arrhythmias occur when the electrical impulses that coordinate the heartbeats do not propagate properly and cause rapid, uncoordinated, weak contractions of the heart 1. Moreover, in the long term, arrhythmias may increase the risk of developing other effects such as stroke and heart failures.

Treatments for arrhythmia can take many different forms, depending on the type and severity of the irregular heartbeat and its cause. Medication treatments for arrhythmia are called chemical cardioversion. The patient receives antiarrhythmic medicine orally or intravenously. But when medications don't offer a solution, an ablation treatment for arrhythmia based on thermal Radio Frequency (RF) can be considered. This ablation technique can be performed on different targets of the heart 2. For example, in atrial/ventricular fibrillation ablation, small scars are intentionally created on the cardiac wall to break up the electrical signals propagation paths that cause the irregular heartbeats. During catheter ablation, a small RF ablator is inserted into the heart, usually through a vein and then small areas of tissue that may be causing or propagating the arrhythmia are necrotized by heating 2. However, RF lesions are not always trans-mural because a good contact between the tissue and the catheter is difficult to achieve 3, and there are also many complications associated with endocardial therapy.

integrates in its middle a 2D US imaging transducers for intraoperative guidance purpose. This imaging device produces an US image perpendicular to the device axis. Classically, the planning of the intervention (ablation path) is defined on a high-resolution preoperative anatomical imaging (CT/MRI). The goal of the cardiologist is then to follow the pre-defined path with the help of the 2D US image alone. But this imaging device used for guidance offers only a very limited view of the anatomy. It seems the important to merge this two information (3D preoperative CT and 2D intraoperative US) by registration in order to improve the guidance of the gesture 7. In a previous work of our team, Sandoval et al. proposed a solution to locate an intraoperative 2D US slice within the 3D CT volume 8. Their method starts from the anatomical hypothesis that the intraoperative 2D US slices are perpendicular to the esophagus axis. So, they reformat the preoperative 3D CT volume into 2D CT slices perpendicular to the esophagus axis. For each 2D CT slice within a candidate zone, they perform a 2D US to 2D CT image-based registration. They use the classical scheme in which an optimizer iteratively searches the transformation (rigid in this case) which gives the highest similarity (Normalized Mutual Information in their case). At the end, the 2D CT slice with the highest similarity after registration is considered as the final result. One of the main drawbacks of this method is that the iterative 2D/2D registration process is too time consuming to be used in the operating room. Recent studies have demonstrated the potential of deep learning methods in directly solving the registration problem without any iterative process. but the main problem of deep learning registration methods is the lack of ground truth for learning. This is why many methods use unsupervised learning 9. However, if the ground truth is known, supervised methods can be considered. Among the class of rigid registration supervised methods, Miao et al. 10 are the first to use deep learning to predict rigid transformation parameters via hierarchical learning. They use a CNN to predict the transformation matrix associated with the rigid registration of 2D/3D X-ray attenuation maps and 2D X-ray images. This approach outperformed the classical image and optimization-based registration approaches in terms of both accuracy and computational efficiency. Recently, Chee et al. 11 use a CNN to predict the transformation parameters used to rigidly register 3D brain MR volumes. In their framework called Affine Image Registration network (AIRNet), the Mean Square Error (MSE) between the predicted and ground truth affine transforms is used to train the network. Salehi et al. 12 use a deep residual regression network, a correction network, and a bivariate geodesic distance-based loss function to rigidly register T1 and T2 weighted 3D fetal brain MRs for atlas construction. However, these methods mainly concern monomodal images.

In this work, we propose a CNN framework for CT/US image registration, using Deep Features representation for a supervised rigid transformation estimation.

METHOD

The heart is a moving organ. However, some characteristics of the cardiac movement allowed us to consider a rigid registration scheme. First, we had at our disposal a Cine CT from a patient's heart composed of 20 volumes at each 5% phase of the RR interval. Second, we are interested in ventricular fibrillation. During its diastolic phase, the ventricle is relatively stationary. The HIFU treatment will be shot in this phase to have a fixed focal point in relation to the organ and thus avoid a dispersion of heat prejudicial to the necrosis of the tissues. So, a quasi-static ventricle pose can be considered. Moreover, on our US system, the acquisition is synchronized with the ECG. Thus, it is relatively easy to create pairs of US/CT images at the same phase and so to consider rigid registration.

In this section, we will introduce the proposed framework for estimating the transformation parameters of a rigid image registration between a preoperative CT slice and an intraoperative US image. In our case and following the approach described in 8, only a 2D rigid transform with three Degrees of Freedom (DOF) -one rotation and 2 translations-has to be estimated.

The main idea is to estimate the registration that best aligns some common characteristics of the images. The information contained in the two images is of a very different nature (gray levels proportional to the X-ray absorption coefficient of the tissues for the CT and information formed by the reflection of waves on surfaces and speckle for US). We must therefore first extract from the two imaging modalities a common information, in our case the shapes of the organs, before performing registration. Therefore, we propose the following process (see Figure 1): (i) Descriptors are extracted from the moving I M and the fixed I F images using Deep Learning. For this, I M and I F are passed through a siamese CNN architecture (ResNet18) consisting of convolutional layers, thus extracting two feature maps f M , and f F which are analogous to dense local descriptors;

(ii) These feature maps are combined in a concatenating layer;

(iii) This corresponding feature maps are the set as input into a convolutional network which directly outputs the parameters set T (two translations, one rotation) of the rigid registration. This framework should be trainable end-to-end for rigid registration task.

Feature extraction

The first step of the framework is feature extraction. Features extraction is a classical tool in deep learning. One common feature extraction technique is to feed the image to a conventional pre-trained neural network and use the representation for that particular image in the intermediate layers of the neural network. We used the ResNet18, which is one of the most efficient standard feature extraction model that can be used in many medical applications 13. This model handles the vanishing or exploding gradient problem when the CNN goes deeper.

Resnet18 can be found implemented in PyTorch. This implementation offers a version with the weights pre-trained for feature extraction on ImageNet, the large benchmark database, Each of our input modalities has its own image characteristics. Thus we passed each of the two images to be registered in its own network. Each of these networks produces a feature map f composed of image descriptors.

Matching

These two feature maps should be combined across images as a single tensor to input it to the rigid transformation parameters estimation network. To achieve this, a concatenation of descriptors along the channel dimensions is performed in the concatenation layer.

Registration network

We will present the network architecture which consists of three blocks of convolutional layers using a kernel size of 5, each followed by batch normalization layers, and a rectified linear unit (ReLU). The last layer is a fully connected layer to estimate the rigid registration parameters. The network expects the concatenated map extracted from moving and fixed images as its input, and directly estimates the parameters (t x , t y , and θ) of the rigid transformation that connects these images. The idea behind this architecture is that the estimation is performed in a bottom-up manner where early convolutional layers vote for candidate transformations, and these are then processed by the later layers to aggregate the votes. The first convolutional layers can also enforce local neighborhood consensus by learning filters which only fire if nearby descriptors in image I M are matched to nearby descriptors in image I F .

Training

We consider a supervised learning scheme, the training datasets include pairwise images of moving CT and fixed US image pairs, and their associated rigid geometric transformation parameters as ground truth GT . The training loss function can be formulated as the L2 norm of the error between the GT (T GT ) and the predicted transformation parameter T Est .

L = α t GT -t Est 2 + β ϑ GT -ϑ Est 2 (1)
With t GT and t Est the translation vector of respectively the ground truth transformation and the estimated one expressed in mm, ϑ GT and ϑ Est the rotation angle expressed in degrees and α and β are weights controlling the balance between the translation and the rotation losses. The choice to use mm for translation and degrees for rotation gave a certain coherence and balance between these parameters. Indeed, an error of 1 degree in rotation leads to a displacement of 1 mm at 60 mm from the center of rotation. Because of this consistency, α and β could be set to 1.

The network is trained by the gradient of the loss function with respect to the estimated rigid parameters (t x , t y , and ϑ). This gradient is then used to minimize the loss function by using backpropagation and Stochastic Gradient Descent.

After training, the network can be applied for registration of unseen image pairs. We implemented the network using PyTorch and we trained it on a NVIDIA TitanX GPU with 10000 iterations.

RESULTS

dataset

Because ground truth cannot be obtained on real data, i.e. we can never ensure that a real US image is perfectly associated to a real CT plane, we decided to simulate US images from CT data.

Our study has been conducted on a public available dataset of twenty contrast enhanced cardiac CT volumes 14,15. All the data were obtained from two state-of-the-art 64-slice CT scanners (Philips Medical Systems, Netherlands) using a standard coronary CT angiography protocol at two sites affiliated to Shanghai Shuguang Hospital. Images were acquired in the axial view, covering the whole heart from the upper abdominal to the aortic arch. The in-plane resolution was about 0.44 × 0.44 mm and the average slice thickness was 0.60 mm.

In these volumes, the esophagus was roughly segmented manually.

From these volumes, we create a set of corresponding pair of CT and US images with known transformation (Figure 2) First we extracted randomly 4000 2D CT oblique cut planes from the 20 CT volumes (200 images per volume). For this we:

1. We choose randomly 4000 initial poses along the esophagus axes within the 20 CT volumes.

2. For each pose we create a new referential by setting some randomly transformations near these initial poses with some translations within ±10 mm and rotations within ±15 degree around each coordinate axis.

3. The x -y plane of this referential will serves as the fixed CT image I CT . The origin of the x -y plane served also as origin of I CT .

4. For each I CT , we simulated the corresponding US images with the method described in 16 that predict the appearance and properties of a B-scan ultrasound image from a probe origin pose, the point spread function of the US device, the acoustical impedance of the tissues and some tissue-adapted distribution of point scatterers. We randomly define the pose of the simulated US origin probe within a range of ±10 mm from the I CT origin and we randomly rotate the probe in a range of ±15 degree around this origin. These two translation and rotation defines the ground truth transformation T GT between the CT and US image.

5. So, for each I CT we get an US image I U S and a transformation ground truth T GT .

From this dataset, the network was trained by selecting the 3600 pairs of corresponding I C T and I U S slices from 18 of the 20 cardiac CT scans. For validation, we used 400 image pairs from the 2 remaining volumes

Evaluation

We compared the registration results obtained by the proposed methods to these obtained by the classical iterative rigid registration method implemented in the SimpleITK Library 17. We used Normalized Mutual Information as similarity measure because it has been found as one of the most suited for our CT/US registration problem 18, We also used 2D Euler transform to presents the spatial mapping of points from the fixed image space to points in the moving image space. We evaluate the performance of the method in terms of computation time and registration accuracy.

Computation time

The mean registration computation time for all the 400 image pairs is now lower than 3 ms for each image pair. which is suitable for a real-time application. For comparison, the classic iterative method takes around 70 seconds to register one image pair.

Transformation estimation error

We compared the parameters of the transformation obtained by our proposed methods with the ground truth (GT ). A transformation is composed by a translation vector t (t x , t y ) and a rotation (vartheta). We evaluate separately the translation errors and the rotation errors between the estimated pose of each of the 400 validation image pairs and their associated GT .

We estimated the translation errors by equation ( 2), where t GT,i and t Est,i are the translation parameters of respectively the GT and the estimated one.

The rotation error has been estimated using the angle difference in degrees estimated between the orientation parameter of the pose of respectively the GT and the estimated rotation angle (equation ( 3)), where ϑ GT,i and ϑ EST,i are the angles that encode the orientation parameter of the pose of respectively the GT and the estimated rotation.

t GT,i -t EST,i 2 (2) |ϑ GT,i -ϑ EST,i | (3)
Figure 3.a shows the boxplots of the 400 translation errors from our CNN-based registration and the classical iterative one. The median translation errors are 1.1 mm using CNN and 1.2 mm using the classical approach.

The boxplots of the rotation errors in Figure 3.b show that the median rotation errors are 2.1 degree using CNN and 2.4 degree when using the iterative classical method.

Target Registration Error (TRE)

The validation can also be done by estimating the registration errors on some fiducial markers. To quantify the error, we defined eight specific feature points (or landmarks) P j in the US fixed images, and we used the two transformations matrices, the estimated T Est and the GT T GT one to project these points in the corresponding CT images: P Est,j = T Est P j and P GT,j = T GT P j . The Euclidean distance between the corresponding projected points P Est,j and P GT,j gives the TRE. Figure 4 shows the boxplot of the TREs for all the 8 fiducial points of all the 400 test images. Quantitative results show a median TRE of 2.2 mm for all the fiducial points of all the 400 test images using CNN, and 2.7 mm using the classical method.

Visual validation

Figure 5 shows a visual comparison between the simulated fixed US image, and the corresponding moving CT image pair. It shows the overlap between the moving CT image and the fixed US image: a) before registration b) after registration with the proposed method. Visually, the results obtained by the proposed method seems to provide a good alignment, this can be seen for example at the probe center, and the bottom of the image on the thoracic chest.

Discussion

From the previous quantitative results, we can conclude that on the one hand, the registration accuracy obtained by CNN is of the same order as that obtained by the classical iterative method. The results obtained by CNN are even slightly better, even if statistically this improvement is not significant. Compared to other methods in the literature, the global target registration error (TRE) of 2.2 mm is on the same range of magnitude as those reported in 8, 10.

On the other hand, CNN allows to strongly accelerate the processing time. The registration between two images takes only 3 ms (instead of 70 s for the classical iterative method). This gain in computation time allows us to consider implementing the 3D CT/2D US registration technique proposed by 8 in clinical practice.

These results were obtained from simulated US images. We are fully aware that there are differences between simulated and real US images (signal attenuation compensation, acoustic shadowing, post processing of real US images...). However, we found in a previous study that a method developed on simulated data performed well on real data 8, 10. We are therefore confident that our method will also work on real data.

CONCLUSION

In this paper, we presented a deep feature learning-based approach for the registration of transesophageal US/CT cardiac images. The results showed a strong improvement in terms of computation time with a promising result in terms of registration accuracy. In future work, we will integrate the features learning approach to a minimallyinvasive HIFU procedure to improve the therapy planning and guidance. We will apply our approach our approach on 2D/3D learning-based registration to refine the estimation of the transesophageal probe pose placement in the 3D preoperative volume. Finally, we will include data from physical phantom and real-patients to evaluate the contribution of our registration scheme to the therapy guidance.

Figure 1 .

 1 Figure 1. The overall of the proposed framework.

Figure 2 .

 2 Figure 2. Datasets creation workflow

Figure 3 .

 3 Figure 3. Box plots of a) the Translation Estimation Errors, b) the Rotation Estimation Errors.

Figure 4 .

 4 Figure 4. Boxplots of the mean Target Registration Errors

Figure 5 .

 5 Figure 5. An example of the registration of an image pair. The overlap between the moving CT image and the fixed US image a) before registration, and b) after registration with the proposed method

ACKNOWLEDGMENTS

This work was part of the CHORUS (ANR-17-CE19-0017) project which have been supported by the French National Research Agency (ANR).


H Labriji
email: hanae.labriji@onera.fr

O Herscovici-Schiller

F Cassaing

Computation of the lateral shift due to atmospheric refraction

Keywords:

Context. Atmospheric refraction causes distortions on the apparent positions of objects in the sky. In complement to the well-known angular offset, we compute here the lateral translation that is also to be considered for nearby objects. Aims. This paper has a dual purpose. We first aim to calculate the lateral shift at each altitude and study its variation according to meteorological conditions and the location of the observation site. We also pay special attention to the chromatism of this lateral shift. And second, we seek to assess the relevance of taking into account Earth's roundness, in relation to the expressions present in the literature in adaptive optics and which have been established neglecting Earth's curvature. Methods. We extract the variation equation of refraction from the geometric tracing of a light ray path. A numerical method along with a model of dry atmosphere are presented, allowing us to numerically integrate the system of coupled equations. In addition to this, we establish three analytic approximations of the lateral shift, one of which is the one already known in the literature. We compare the three approximations to the numerical solution. Results. Regarding the numerical computation of the lateral shift, we studied its evolution linked to the variation of temperature and pressure at the observer level. Concerning our second order approximation, we found it more precise (accuracy better than 1% for zenith angles up to 75 • ) and less pessimistic than previous estimates. We establish that the first order approximation is largely sufficient to estimate the impact of the lateral shift on adaptive optics systems. Finally, we found that the lateral shift must be considered when compensating atmospheric refraction for the observation of nearby objects such as meteors space debris.

Introduction 1

Derivation of the lateral shift expressions 2 2.1 Differential equations from geometric analysis . . 2 2.2 Numerical integration of the refraction integral . 4 2.3 Second order approximation of the local equations 5 2.4 Analysis of the terms involved in the second order approximation . . . . . . . . . . . . . . . . . 6

Accuracy of the second order approximation and evaluation of the lateral shift variability A A model for the dry atmosphere 13

Introduction

Due to the decrease of air refraction index according to altitude, light rays entering the atmosphere or emitted from Earth are bent, respectively, towards the ground or towards the zenith. One consequence of this phenomenon, illustrated in Fig. 1, is the difference between the apparent direction of the stars (position number (1)) and their true directions (position number (0)). This deviation is called the angle of refraction R a and depends on the observing wavelength, which causes a visible stretching of the stars also known as the blur effect. Atmospheric refraction was very early recorded in the literature, Aristotle mentions the vertical stretching at the horizon in his book Meteorology [START_REF] Aristotle | Works. Translated into English under the editorship[END_REF]) while the first instrumental proof was made by Johannes Schöner and reported in a monograph written by [START_REF] Regiomontanus | Scripta de Torqueto, Astrolabio armillari, Regula magna Ptolemaica[END_REF]. The famous Danish astronomer Tycho Brahe was one of the first persons to measure the effect of refraction using the apparent position of sun at summer and winter Earth's solstices (Mahan 1962). Later, the angle R a has been the subject of several studies, authors such as Radau (1882) and [START_REF] Saastamoinen | [END_REF] have approached the refraction angle with a series of approximations while others integrated numerically the variation equations of the light path [START_REF] Hohenkerk | The Computation of an Angular Atmospheric Refraction at Large Zenith Angles[END_REF][START_REF] Seidelmann | [END_REF]Stone 1996;Wittmann 1997;[START_REF] Kristensen | [END_REF][START_REF] Auer | [END_REF]Nauenberg 2017). There are also several equation-based refraction tables, such as the one made in the observatory of Pulkovo (Kurzyńska 1988).

The purely angular correction R a is well adapted to astrometry, i.e. when compensating the position of an object optically located at infinity such as a star. But there are a number of cases where it is necessary to take into account the shift b of the ray with respected to its corrected path, pictured in Fig. 1.

By way of illustration, the shift compensation is relevant in the three following examples. A first case is estimating the position of an object near Earth using its apparent position in the sky, this is typically used in meteor tracking photography networks (McCrosky & Posen 1968;Jeanne et al. 2019;Gardiol et al. 2021;[START_REF] Colas | Proceedings of the International Meteor Conference[END_REF]Colas et al. , 2020)). Since the object is generally observed at an altitude larger than 20 km, its apparent position is affected by atmospheric refraction. In meteorite tracking networks, refraction is often included in the dispersion correction of cameras by taking as reference stars located at infinity to calibrate the all-sky cameras (Borovicka et al. 1995). Yet, the correction to be applied is more than just R a and has been highlighted by McCrosky & Posen (1968) as the "refractive parallax". To illustrate, we consider again Fig. 1 where we draw the path of a light ray emitted by a meteorite when it enters Earth's atmosphere. The apparent position of the meteorite corresponds to the apparent position of the faraway star. When this direction is corrected by only the refraction angle R a , there is clearly an error equal to the shift b in the estimated position of the observed object.

A second example of the shift compensation usefulness is the observation of Earth from space using remote sensing. This time, the projection of the lateral shift on the ground is used to compensate refraction when estimating the position of the generated data (Yan et al. 2016;Li et al. 2016). This distance is denoted by ∆ in Fig. 1. Also, since atmospheric wavefront distortions are much less critical when observing Earth from space than space form Earth, Earth observation satellites usually operate at large zenith angle and thus, suffer more from the effect of refraction.

The third case is adaptive optics, where the optical path difference according to the observation wavelength will induce an error in the achromatic correction of atmospheric turbulence. This is a form of "chromatic anisoplanatism" [START_REF] Devaney | [END_REF]Sasiela 1992) sometimes called the "chromatic shear" (Nakajima 2006). For instance, adaptive optic systems correct atmospheric turbulence using achromatic devices, except that light rays of different colors will pass through different paths in the atmosphere and thus undergo different phase distortions. Although the angular dispersion is usually compensated fully by an atmospheric dispersion corrector, the impact of refraction on atmospheric turbulence is, to date, mostly not taken into account (as pointed by van den Born & Jellema (2020)).

The calculation of the lateral shift b has been identified as a difficult problem by McCrosky & Posen (1968) and there have been a few early attempts to approximate it [START_REF] Schmid | The Influence of Atmospheric Refraction on Directions Measured to and from a Satellite[END_REF]. Within the framework of adaptive optics, Wallner has written several papers on this issue. The most relevant here are [START_REF] Wallner | Imaging Through the Atmosphere[END_REF] and [START_REF] Wallner | [END_REF]). In [START_REF] Wallner | Imaging Through the Atmosphere[END_REF], one can find an expression of the lateral shift whose derivation assumes that the Earth is flat:

b(h) ∝ tan z ∞ cos z ∞ P(h) gρ S , (1)
where z ∞ refers to the zenith angle of the unrefracted ray, P is the atmospheric pressure at altitude h, g is the standard gravity and ρ S is a reference value for air density.

Based on Eq. 1, some authors have continued his work to evaluate losses for extremely large telescopes (ELTs) [START_REF] Owner-Petersen | [END_REF][START_REF] Owner-Petersen | Advances in Adaptive Optics II[END_REF][START_REF] Jolissaint | 1st AO4ELT conference -Adaptive Optics for Extremely Large Telescopes[END_REF]. Several other papers seem to have independently found and used the same expression of the lateral shift (Sasiela 1992;Nakajima 2006;[START_REF] Devaney | [END_REF].

Our work is dedicated to the derivation of the lateral shift in the spherical Earth model. We begin by rigorously writing the differential equations from which we deduced three approximations of the lateral shift, providing the means to calculate them using only the meteorological conditions of the observation site and the distribution of air in the atmosphere. In Sect. 3, we study the accuracy of each approximation as a function of the zenith angle. Then, the system of coupled equations is solved for a set of various parameters in order to investigate the dependence of the lateral shift on observing conditions. Finally, we study in Sect. 4 two applications of the lateral shift. We start by studying in Sect. 4.1 the effect of numerical integration on the "chromatic shear" values in adaptive optics, and then we focus on the effect of lateral shift for the observation of nearby objects (Sect. 4.2).

Derivation of the lateral shift expressions

Differential equations from geometric analysis

We suppose a star at infinity, we take a ray that reaches the observer at point Ω and trace back its optical path in a spherical model of Earth's atmosphere as depicted in Fig. 2. The light ray passes through point M at the altitude h and point N at h + dh. The hypothetical path in absence of atmosphere is represented by a dotted line. A is the orthogonal projection of M on the unrefracted ray, while A is the orthogonal projection of N on (MA ).

We note:

-R T Earth's radius, h the altitude of the current point, z the angle that the ray makes with the direction of the zenith at the observer, -ζ the local zenith angle at each altitude, s the length along the optical path, n the local refractive index.

The subscript 0 always refers to a quantity measured at ground level (i.e. at the observer level) while the subscript ∞ refers to a quantity at infinity. z 0 is the apparent zenith angle at the observer level and z ∞ is the angle between the unrefracted ray and the zenith at Ω. The classical refraction angle is denoted R a and: The variation of the local zenith angle at each altitude h satisfies the Bouguer invariant and depends on the local refractive index as follow:

R a = z ∞ -z (2)
R T n 0 sin ζ 0 = (R T + h) n sin ζ.
(3)

This formula stems from Snell Descartes' law, its proof can be found in chapter 6 of [START_REF] Kovalevsky | Fundamentals of Astrometry[END_REF].

In the following we derive the equations of variations of the main variables describing the optical path of the light ray: θ, ζ and h. Some of those are already known in the usual computation of the angle of refraction z 0z ∞ [START_REF] Hohenkerk | The Computation of an Angular Atmospheric Refraction at Large Zenith Angles[END_REF]. The goal here is to express especially the lateral shift b(h) according to given parameters such as the zenith angle. As already highlighted in many references [START_REF] Kovalevsky | Fundamentals of Astrometry[END_REF], the main difficulty of this calculation in spherical atmosphere is the continuous change of the zenith direction (ζ) along the path of the ray. Using the variable z, which represents the angle that the ray makes with the zenith at the point of observation, allows us to overcome this difficulty.

First, differentiating the Bouguer formula, one gets:

tan ζ dh R T + h + dζ = -tan ζ dn n , (4)
Then in Fig. 2, one has in the triangle OK M:

z = θ + ζ. (5)
Moreover, in the triangle MHN (where MH is a circle arc):

tan ζ = (R T + h) dθ dh , (6)
which implies:

dz = dζ + dh R T + h tan ζ. (7)
which, with Eq. 4, gives:

dz = - dn n tan ζ. (8)
To obtain numerically the values of z, θ and ζ along the optical path, it is necessary to integrate the variation equations according to a chosen integration variable. There are several possible integration variables. Knowing n and dn/dh, the choice of h as the integration variable is the most intuitive, but reveals a singularity near the horizon because of tan (ζ) in Eq. 8 (the tangent function is not defined in π/2). [START_REF] Auer | [END_REF] recommend using z as the integration variable, but this adds a computational step and makes a singularity appear at the zenith [START_REF] Hohenkerk | The Computation of an Angular Atmospheric Refraction at Large Zenith Angles[END_REF]. We choose here to use s as recommended by van der Werf (2008) because it is the one that avoids singularities at both ends.

Considering again the triangle MHN, we can connect the infinitesimal variations of s to h using:

dh = cos (ζ) dh. ( 9
)
From Eqs. 8 and 9, we got:

dz = - sin (ζ) n dn dh ds, (10)
and from Eqs. 6 and 9, we got also:

dθ = sin (ζ) R T + h ds. ( 11
)
As a consequence, we merge the Eqs. 10 and 11 with Eq. 5 to obtain:

dζ = - sin (ζ) n dn dh . ( 12
)
We deal next with the lateral shift b. Since A is the orthogonal projection of N on (MA ), one can write:

b(h) = b(h + dh) + AM. ( 13
)
Then, considering the triangles AMN et BNK, one has:

AM = sin (z ∞ -z) ds = -db, (14)
and thus:

db = sin (z -z ∞ ) ds (15)
Finally, combining Eqs. [9-12] and Eq. 15, we end up with a system of five first order coupled non linear ordinary differential equations:

dh ds = cos (ζ), ( 16a
) dz ds = - sin (ζ) n dn dh , ( 16b
) dθ ds = sin (ζ) R T + h , ( 16c
) dζ ds = - sin (ζ) n dn dh - sin (ζ) R T + h , ( 16d
) db ds = sin (z -z ∞ ). ( 16e
)
Of the five variables, one can choose to omit one of z, θ et ζ because they are linked by Eq. 5.

Numerical integration of the refraction integral

In order to solve the system of Eqs. 16, it is necessary to provide n and dn/dh over the entire height of the atmosphere. Since the true values of the refractive index are not easily obtained, it is common to choose a particular model of the atmosphere. We make the assumption that the atmosphere is a perfect gas of low density and is in hydrostatic equilibrium. It follows the Gladstone Dale relation (Gladstone & Dale 1863):

n -1 = κρ. ( 17
)
where κ is a constant which depends on the observation wavelength λ, and ρ is the local density of air. In order to express the density ρ as a function of altitude, it is sufficient to get a temperature profile. We choose to consider a constant temperature gradient in the tropophere:

ω = dT /dh = -6.5 K/km (18)
and a constant temperature beyond the tropopause (i.e. the temperature gradient is nil). This will allow us to adapt the evolution of the temperature according to its value T 0 measured at ground level. Thus, our inputs are the temperature and pressure at ground level (respectively T 0 and P 0 ) and the relative humidity (H R ) considered as constant in the troposphere and nil above it. This model truncates the US76 standard atmospheric profile that is based on a piecewise linear temperature for nine atmospheric layers (COESA 1976). The choice to simplify it and reduce it to only two layers is motivated by the fact that most of the atmospheric refraction occurs in the troposphere, since the air beyond it is very sparse. Many authors have made this approximation for the computation of the angle of refraction [START_REF] Hohenkerk | The Computation of an Angular Atmospheric Refraction at Large Zenith Angles[END_REF]Néda & Volkan 2002;[START_REF] Auer | [END_REF]. Details on the dry atmosphere model are given in Appendix A. Now that we have chosen a model for the atmosphere and that we are therefore able to compute the values of n and dn/dh at each altitude, we present next the numerical arguments of solving the system of coupled ordinary differential Eqs. 16.

First, we notice that db/ds variation equation also requires knowledge of z ∞ . Therefore, it is necessary to proceed in two separate steps:

1. solving the first equation system ) to calculate the refraction angle R a , 2. and then solving the complete system (Eqs. [16a -16e]) to derive the lateral shift b.

Both computations are made using a fourth order Runge-Kutta method, and in the reverse direction of the optical propagation. We place ourselves in the point of view of the observer and take the reference of the lateral shift at the ground. This way, we consider the following initial state:

                 h 0 z 0 0 z 0 0                  =                  h z θ ζ b                  h=h 0 . ( 19
)
10 0 10 1 10 2 10 3 Integration step (m)

10 -9 10 -8 10 -7 10 -6 10 -5 10 -4 10 -3 ∆R a /R a z 0 = 65 • z 0 = 75 • z 0 = 85 • Fig. 3. ∆R a /
R a versus the integration step.

We relate the integration limit in h and s using the approximation :

S max = H max cos z 0 . ( 20
)
This expression is not appropriate close to the horizon since the refraction angle R a grows and the approximation of 1/cos ζ requires higher orders. However, this value of S max is an upper bound for the effective path length of the ray in the atmosphere of width H max .

Then, in order to set an integration step ds, we compute the angle of refraction for several steps starting from ds = 2000 m. This step is divided by two at each iteration, until the relative variation between two successive solutions becomes less than 10 -10 . A very low value of the minimum relative deviation was consciously chosen in order to obtain the most reliable reference for evaluating the error made by the numerical integration. It goes without saying that the other inherent errors in the atmospheric model are much larger and will prevent such accuracies from being achieved.

The relative deviation between the refraction angle at each ds and its value for the smallest integration step ds min is plotted in Fig. 3:

∆R a /R a = R a (ds) -R a (ds min ) R a (ds min ) , (21)
and at three different zenith angles 65 • , 75 • and 85 • . We notice in Fig. 3 that the relative error decreases with ds. Also, the decay slope is lower for z 0 = 85 • , but in all cases the relative error reaches values less than 0.01% for all the integration steps lower than 100 m. In the following, we choose ds = 100 m. Also, we solve each system separately in the troposphere and in the stratosphere because of the discontinuity of the refractive index derivative at the tropopause (due to the discontinuity of temperature gradient). To do so, we make a change in the integration step at the interface. This way, the last integration step in the troposphere is the only one to be changed.

Computer programming is made in Python, the source code (has been/will be) shared on the collaborative platform GitHub at the following web address to Add.

Although the numerical integration of the differential equations can be done quickly using a modern computer, an approximate analytical formula can often be handy. That is the purpose of Sect. 2.3.

While we previously chose s the distance along the optical path as the integration variable, we have to define t he integration limit as a function of the maximum height of Earth's atmosphere H max : the altitude at which the density of air can be regarded as insignificant. The value of H max is often chosen equal to 80 km in the literature [START_REF] Hohenkerk | The Computation of an Angular Atmospheric Refraction at Large Zenith Angles[END_REF]van der Werf 2003).

Second order approximation of the local equations

The system we have achieved in Sect. 2.1 can not be solved analytically, nor does it allow to reflect the influence of the main parameters (Earth's roundness, density distribution, temperature and pressure at the observer) on the value of the lateral shift. In order to clarify this, we approximate the value of the lateral shift b in a similar way to Laplace's formula for the refraction angle (Radau 1882). Laplace's formula refers to an expression that allows the refraction angle to be calculated from the apparent zenith angle and the weather conditions at ground level only, it states that:

R a = α (1 -β) tan z 0 -α β - α 2 tan 3 z 0 ( 22
)
where α is defined by:

α = A D (λ) P 0 T 0 , ( 23
)
where P 0 and T 0 are respectively temperature and pressure at the observation site, and A D is the refractivity of dry air defined in Eq. A.3 (α still depends on the observation wavelength). And

β = P 0 /(h 0 + R T )ρ 0 .
To begin with, we take advantage of the following:

h R T , (24) n -1 1, and (25) sin (z -z ∞ ) z -z ∞ . (26)
The first two approximations are perfectly justified, on the one hand h/R T has as an upper bound H max /R pol ≈ 10 -2 (where R pol = 6356.75 km is the polar radius of Earth, and H max is the height of the atmosphere defined in Eq. 20) and on the other hand the value of refraction index has its maximal value at ground level and does not exceed 10 -3 + 1. The sine approximation is quite reasonable too because the angle of refraction is always very small and does not exceed 30 arc-minutes, even when considering real weather data as in the study carried by Nauenberg (2017).

Using the parameter α defined in Eq. 23, the refraction index is written:

n = 1 + α T 0 P 0 P T . (27)
And since atmospheric air is assumed to behave according to the ideal gas law, one also has:

n = 1 + α ρ, ( 28
)
where we denote by ρ the normalized density ρ/ρ 0 . Knowing Eqs. 23 and 25, we have α 1 as well.

The approximation of the lateral shift will be written as a function of the main characteristic lengths of the atmosphere defined below:

L 1 (h) = ∞ h ρ dx, ( 29a
) L 2 (h) = ∞ h ρ2 dx, ( 29b
) L b (h) = ∞ h x ρ dx /L 1 (h). ( 29c
)
Using the hydrostatic equilibrium equation, L 1 can also be written as:

L 1 (h) = P(h) gρ 0 . (30)
Furthermore, L b (h) is a length that can be interpreted as the altitude of the barycenter of an air column.

In order to get the approximation, we start by integrating Eq. 15 between h 0 (the altitude of the observation site) and the limit of Earth's atmosphere:

b(h = h 0 ) -b(h → ∞) = h 0 ∞ sin (z -z ∞ ) ds. ( 31
)
The reference of b can be chosen either at ground level or at infinity, depending on the purpose of the calculation. We choose to make the following calculations with b(h → ∞) = 0.

In order to use the model of the atmosphere exposed in the Appendix A, we perform a change on the integration variable:

b(h = h 0 ) = h 0 ∞ sin (z -z ∞ ) cos ζ dh. ( 32
)
And we restrict ourselves to the directions of observation well above the horizon to avoid divergences. Then, we perform the integration of Eq. 8 between an altitude h and ∞ along with a Taylor expansion of the integrand for small α and small h/R T following [START_REF] Kovalevsky | Fundamentals of Astrometry[END_REF]. And so we get the following approximation of zz ∞ :

z(h) -z ∞ = -tan z 0 α ρ -α 2 ρ2 + α 2 ρ -hα ρ/R T -αL 1 (h)/R T + tan 3 z 0 α 2 ρ2 /2 -α 2 ρ + hα ρ/R T + αL 1 (h)/R T . (33)
Then we expand the cosine term in Eq. 32 using the Bouguer's invariant (Eq. 3):

1 cos ζ = 1 1 -sin 2 ζ , (34a)
= (h + R T ) n/(n 0 R T ) [(h + R T ) n/(n 0 R T )] 2 -sin 2 z 0 . (34b)
In addition to this, since α 1 and h R T , one has:

(h + R T ) n R T n 0 = 1 + h R T 1 + α ρ 1 + α , (35a)
= 1 + h R T (1 + α ρ) (1 -α + o(α)) , (35b)
= 1 + h R T + α( ρ -1) + o(α + h R T ). ( 35c
)
And based on this, we develop Eq. 34b to the second order with respect to (h + R T ) n/R T n 0 . One gets:

1 cos ζ = 1 cos z 0 - tan 2 z 0 cos z 0 (h + R T ) n R T n 0 -1 + 3 2 1 cos 5 z 0 - 1 cos 3 z 0 (h + R T ) n R T n 0 -1 2 . ( 36
)
As mentioned before, we note that this approximation is less accurate when the zenith angle z 0 approaches π/2 (i.e. when looking close to the horizon), the reason is that 1/ cos z 0 is no longer defined when z 0 = π/2.

We perform then the product of the two approximations: Eqs. 33 and 36. This gives a second-order approximation b (2) of the lateral shift in a dry atmosphere:

b (2) = tan z 0 cos z 0 ∞ h 0 α ρ -α 2 ρ2 + α 2 ρ -αh ρ/R T -αL 1 (h)/R T dh - tan 3 z 0 cos z 0 ∞ h 0 3α 2 ρ2 /2 -2α 2 ρ + 2αh ρ/R T + αL 1 (h)/R T dh. ( 37
)
Using an integration by part, one can also prove that:

∞ h 0 L 1 (h) dh = - h 0 P 0 gρ 0 + ∞ h 0 h ρ dh. ( 38
)
We then get a simple approximation of the lateral shift depending on the moments of the atmosphere defined in Eq. 29a :

b (2) = A sh tan z 0 cos z 0 -B sh tan 3 z 0 cos z 0 , (39)
where:

A sh = αL 1 (h 0 )+ αh 0 R T L 1 (h 0 )-2αL b (h 0 ) L 1 (h) R T -α 2 L 2 (h 0 ) -L 1 (h 0 ) , (40)
and

B sh = - αh 0 R T L 1 (h 0 )+3αL b (h 0 ) L 1 (h 0 ) R T +α 2 3L 2 (h 0 )/2 -2L 1 (h 0 ) . (41)
Besides, using the atmosphere model exposed in Appendix A we can explicitly express the moments L 2 (h 0 ) and L b (h 0 ) as a function of T 0 , H t , and the ratio of T t and T 0 :

r := T t T 0 . ( 42
)
We get:

L 2 (h 0 ) = γL 1 (h 0 ) 2γ -1 - L 1 (h 0 ) 2(2γ -1) r 2γ-1 , and (43)
L b (h 0 ) = h 0 1 - r γ γ + 1 + γ γ + 1 L 1 (h 0 ) 1 -r γ + γ + 1 γ r γ+1 + r γ H t γ + 1 . ( 44
)
The evolution of L 1 L 2 and L b as a function of the altitude at the observation point h 0 is shown in Fig. 4. In order to obtain this figure, we have changed both the altitude and the temperature and pressure values at the observation site. The evolution of P(h 0 ) and T (h 0 ) is based on the model of the dry atmosphere exposed in Appendix A using the standard values of P(h 0 = 0) and T (h 0 = 0):

T 0 = 273.15 K, P 0 = 1000 hPa, (45)
L 1 (h 0 ), L 2 (h 0 ) and L b (h 0 ) (m) L 1 L 2 L b

Analysis of the terms involved in the second order approximation

In the approximate analytic formula for the lateral shift established in Sect. 2.3, one can distinguish three different terms:

a term proportional to α, another term proportional α/R T which is due to Earth's roundness, one last term proportional to α 2 .

We call these three distinct components respectively b

(1) , b (2)
R and b (2) 2 :

b (1) (h 0 ) = α tan z 0 cos z 0 L 1 (h 0 ), (46a) b (2) R (h 0 ) = α R T tan z 0 cos z 0 L 1 (h 0 ) h 0 -2L b (h 0 ) + α R T tan 3 z 0 cos z 0 L 1 (h 0 ) h 0 -3L b (h 0 ) , (46b) b (2) 2 (h 0 ) = α 2 tan z 0 cos z 0 L 1 (h 0 ) -L 2 (h 0 ) -α 2 tan 3 z 0 cos z 0 3 2 L 2 (h 0 ) -2L 1 (h 0 ) . (46c)
b (1) is the first order approximation of the lateral shift. It is also the expression of the lateral shift as given in many references [START_REF] Wallner | Imaging Through the Atmosphere[END_REF]Sasiela 1992;Nakajima 2006;[START_REF] Devaney | [END_REF] where the calculations was performed under the assumption of a flat Earth and a parallel atmosphere.

In order to capture the influence of each term, we plot on Fig. 5 the evolution in absolute values of each of the three terms according to zenith angle. For the calculation of the multiplying factor alpha, we choose the wavelength λ 0 = 550nm as a reference value, it is the middle of the visual spectral band bande V c'est bien pour visible ?.

In Fig. 5, we observe that the main contribution comes from the flat Earth term and is followed by Earth's roundness. The dominance of the term due to Earth's roundness drives us to define a third estimator of the lateral shift denoted b (3/2) , and composed of the two terms b (1) and b (2) R such that:

b (3/2) = b (1) + b (2) R . (47)
The final expression in Eq. 39 efficiently decouples the contribution of the atmosphere to the lateral shift from the contribution of the apparent zenith angle and the observation wavelength. It also spares the user to have to integrate an integrand varying with the zenith angle. Thus, we have somehow established a Laplace formula for the lateral shift. Nevertheless, there is a major difference between the two formulas because of the two moments L 2 and L b in the expression of b (2) . Indeed, these two lengths depend on the composition of the atmosphere and the distribution of air within it. In the following section, we describe in detail the terms that are part of our approximation. This third approximation, whose order is between the first and second orders, has the advantage of decoupling the wavelength from the zenith angle. For instance, the lateral shift is separable into functions of λ and z 0 :

b (3/2) (z 0 ) =α tan z 0 cos z 0 L 1 (h 0 ) + L 1 (h 0 ) R T h 0 -2L b (h 0 ) + L 1 (h 0 ) R T tan 3 z 0 cos z 0 h 0 -3L b (h 0 ) , (48)
α being a function of the observation wavelength. We plot on Fig. 6 the value of:

∆b (3/2) (z 0 , λ) b (3/2) (z 0 , λ 0 ) = b (3/2) (z 0 , λ) -b (3/2) (z 0 , λ 0 ) b (3/2) (z 0 , λ 0 ) (49)
as a function of λ. This quantity depends only on the wavelength λ, and is in fact simply equal to:

α(λ) -α(λ 0 ) α(λ 0 ) . ( 50
)
This way, one can estimate the value of the lateral shift for any wavelength λ and any zenith angle z 0 using this curve and the value of the lateral shift at λ 0 and the right zenith angle z 0 .

Accuracy of the second order approximation and evaluation of the lateral shift variability

Approximation error

In order to quantify the accuracy of the previously established approximations, we calculate the numerical value of the lateral shift resulting from the integration of Eq. 16 for various zenith angles and compare these values with the three estimators b (1) , b (3/2) and b (2) which also depend on z 0 , all at the standard conditions of temperature and pressure (SCTP), at sea level and with an observation wavelength equal to λ 0 . We plot the absolute and relative errors of the three approximations of the lateral shift respectively in Figs. 7 and8:

(X) a (z 0 ) = b (X) (h 0 = 0, z 0 ) -b 0 (z 0 ) , ( 51
) (X) r (z 0 ) = b (X) (h 0 = 0, z 0 ) -b 0 (z 0 ) b 0 (z 0 ) , ( 52
)
where X stands respectively for 1, 3/2 and 2. It is clear that the three approximations deviate from the numerical solution as the zenith angle increases. This is consistent with the approximation that has been made in Eq. 36 and that is also implicitly present in Laplace's formula and therefore in Eq. 33. From Fig. 8, we can see that the second order approximation has less than 1% relative error up to a zenith angle of 75 • while the first-order approximation satisfies this requirement only up to z 0 = 55 • . Concerning the approximation b (3/2) , it is right in between the two "complete" order approximations, it suffers 10 times less relative error than order 1 and thus achieves less than 1% relative error up to z 0 = 70 • .

Influence of temperature and pressure

This part focuses on the influence of weather conditions on the value of the lateral shift. Thus, we start by plotting in Fig. 9 the evolution of the lateral shift at ground level as a function of the apparent zenith angle using the numerical method described in Sect. 2.2. The results were computed for the standard conditions of temperature and pressure (SCTP):

T 0 = 273.15 K, P 0 = 1000 hPa, ( 53
)
and without humidity: H R = 0%. We can see that the lateral shift increases with the apparent zenith angle, it tends towards 0 10. Lateral shift versus ground temperature for various zenith angles at λ 0 and using the numerical integration of Eqs. 16. 11. Lateral shift versus ground pressure for various zenith angles at λ 0 and using the numerical integration of Eqs. 16.

z 0 = 5 • z 0 = 25 • z 0 = 45 • z 0 = 65 • z 0 = 75 • Fig.
z 0 = 5 • z 0 = 25 • z 0 = 45 • z 0 = 65 • z 0 = 75 • Fig.
since each light ray at wavelength λ is bent by a different refraction angle R a (λ). This issue is more pronounced when looking towards the horizon. In this section we study the evolution of the lateral shift as a function of wavelength and with respect to λ 0 = 550 nm. We plot on Fig. 12 the lateral shift at h 0 = 0 as a function of wavelength and for several apparent zenith angles.

∆b 0 (z 0 , λ) = b 0 (z 0 , λ) -b 0 (z 0 , λ 0 ). (54)
The results are computed for the standard conditions of temperature and pressure, without humidity, and over a range of wavelengths from 0.4 µm to 4.2 µm. Since Ciddor's refractivity formula (Ciddor 1996) is restricted to a maximum wavelength of 1.49 µm, we have covered the remaining interval using Mathar (2007) formulas. We observe that the lateral shift decreases with the wavelength for any zenith angle. In order to quantify the relative extent of this decrease, we plot in Fig. 13 the relative difference in b(λ) with respect to its value at λ 0 and for an apparent zenith angle z 0 = 65 • , 75 • and 85 • :

∆b 0 (z 0 , λ) b 0 = b 0 (z 0 , λ) -b 0 (z 0 , λ 0 ) b 0 (z 0 , λ 0 ) . ( 55
)
We see that the variation of the lateral shift depending on the observation wavelength reaches at most 3%.

Fig. 9. Lateral shift at SCTP versus apparent zenith angle at λ 0 = 550nm and computed using the numerical integration of the lateral shift equations 16.

when z 0 tends towards 0 and increases with a very steep slope near z 0 = 90 • , which matches perfectly with intuition and the evolution of the refraction angle R a . Indeed, the greater the zenith angle, the greater the angle of refraction and the greater the ray path through the atmosphere. We note also that the lateral shift reaches significant values beyond 2 0 • (b ≥ 1 m), and can exceed tens of meters beyond 60 • . Then we take a closer look at the distinct effects of temperature and pressure. We plot on Fig. 10 the evolution of b as a function of T 0 for several apparent zenith angles, and on Fig. 11 the influence o f P 0 . We can see on Fig. 10 that the lateral shift decreases on a small scale with temperature but does not exceed 1% of its value at T 0 = 253.15K. In contrast, pressure has the opposite influence, indeed the lateral shift increases slightly with the pressure. In the defined pressure range, the increase of the pressure does not exceed 7% of its value at P 0 = 98000 Pa, regardless of the zenith angle. The greater the apparent zenith angle, the more visible both effects are. These two variabilities are in line with our expectations, since air density increases with pressure and decreases with temperature.

Dependence on observation wavelength

Since the refractive index depends on the observation wavelength λ, the effects of refraction and more specifically the lateral shift also depend on λ. This results in a chromatic blur effect (illustrated in Fig. 1) that appears when observing a celestial body We can also notice in Fig. 13 that the curves of ∆b 0 /b 0 are almost identical for the three zenith angles considered, it is resulting from the separability of the approximation b (3/2) of the lateral shift established in the paragraph 2.4.

z 0 = 5 • z 0 = 25 • z 0 = 45 • z 0 = 65 • z 0 = 75 •

Application to meteor trajectography and adaptive optics

4.1. Relevance of the second order approximation of the lateral shift for adaptive optics This section focuses on the impact of the lateral shift in the performance of adaptive optic systems. Indeed, as the turbulence corrector is achromatic, there is a chromatic residual phase that is not corrected and is due to the fact that light rays of different colors go through different layers of turbulence as well. We consider an adaptive optics system where the measurement and analysis of the wavefront is done at a wavelength λ WFS = 550 nm, while the scientific imaging after compensating atmospheric turbulence is done at a wavelength λ SCI freely selected. As indicated in [START_REF] Wallner | Imaging Through the Atmosphere[END_REF] and Nakajima (2006), the values of the lateral shift at the two wavelengths (λ SCI and λ WFS ) all along the ray path allow us to compute the wavefront phase variance associated to the chromatic shear, by the use of a specific and known turbulence profile. Fig. 14. Difference between the lateral shift at λ WFS and various λ SCI versus zenith angle, computed using the numerical integration of Eq. 16.

|∆b 0 | (m) B R I H K L M N
Our goal here is not to duplicate the mathematical developments present in the previously cited papers, but only to discuss the relevance of calculating the lateral displacement of the ray more precisely, i.e. using numerical computation instead of the first order approximation.

To begin with, we plot on Fig. 14 the difference between the value of the lateral shift at λ WFS and its value at λ SCI as a function of the observation zenith angle. Several values of λ SCI are used, these are the middle wavelengths of the spectral bands B (445 nm), R (658 nm), I (806 nm), H (1630 nm), K (2190 nm), L (3450 nm), M (4750 nm), and N (10500 nm). We define ∆b 0 as:

∆b 0 = b 0 (z 0 , λ WFS ) -b 0 (z 0 , λ SCI ) (56)
The lateral shift was on purpose calculated over the entire thickness of the atmosphere, so as to simplify the computation and not consider any specific turbulence profile. We observe on this figure that the deviation ∆b 0 increases both as a function of the zenith angle and the wavelength. The greater the difference between λ WFS and λ SCI , the greater is ∆b 0 , because of the monotonic nature of the lateral shift as a function of the wavelength. We also notice the stationarity that we see in Fig. 12, the curves of ∆b 0 with a λ SCI greater than 1630 nm are almost identical. We note that the B band is located between the R and I spectral bands due to the absolute value.

Next, we take a closer look at the effect of numerical integration. According to equation 10 of Nakajima (2006), the phase variance of the wavefront is proportional to (∆b 0 ) 5/3 . We plot on Fig. 15 the ratio between ∆b 0 , calculated using the numerical integration of the system of Eqs. 16, and ∆b (1) 0 , calculated using the first-order approximation (Eq. 46a), all to the power 5/3 and for λ SCI = 1630 nm. We notice on Fig. 15 that the first order approximation overestimates the values of the lateral displacement of the rays as the ratio is less than 1 regardless of the zenith angle. Also, the ratio studied decreases when the angle z 0 increases, but is still very close to 1, reaching 99.1% at z 0 = 50 • .

Besides, we may wonder about the significance of this decrease in terms of sky coverage. To find out, we have calculated the zenith angles (z 1,i 0 ) for which the values of ∆b 0 (computed solving Eq. 16) is equal to the ones found with the first order approximation for a sequence of zenith angles (z i 0 ). The difference between the two sequences of angles is plotted on Fig. 16 as a function of the initial zenith angles (z i 0 ), this ing position of the object from pictures taken at different locations and by several cameras. As previously described in the introduction, the correction by the refraction angle which is intrinsic in the change of the camera's reference frame, is not enough to compensate for atmospheric refraction and there is then an error exactly equal to the lateral shift. If we want to compare the impact of these two corrections, the one by the refraction angle and the one by the lateral shift, we can write the compensation by the lateral shift as a correction through the angle σ defined in Fig. 1 and with:

z 0 = 5 • z 0 = 25 • z 0 = 45 • z 0 = 65 • z 0 = 75 • z 0 = 80 • z 0 = 85 •
sin σ = b(H) D . ( 57
)
where D is the true distance between the object and the observer and H is the altitude of the object.

In order to know whether this additional correction is significant or not in relation to the main refraction angle, we plot on Fig. 17 the evolution of the ratio σ/R a according to altitude. Equation 58 gives the distance D as a function of the object's altitude H, Earth's radius R T and the angle θ defined in Sect. 2.2.

D 2 = R 2 T + (R T + H) 2 -2R T (R T + H) cos θ. (58)
Based on the reported observations of meteors (Jeanne et al. 2019;Gardiol et al. 2021), 80 km is the typical altitude at which a meteor begins to be observed, and it disappears around an altitude of 20 km as it comes apart in the atmosphere. The observation site is assumed at see level, under the standard conditions for temperature and pressure and with an observation wavelength λ 0 .

We notice on Fig. 17 that the closer the object gets to Earth and the more important the correction of angle σ becomes, it even reaches 1/3 of the correction angle R a at H = 20 km. This leads to the fact that the correction of the lateral shift must be taken into account to improve the performance of the detection systems. Also, we see that the ratio depends only slightly on the zenith angle for small to moderate angles and then decreases strongly for large zenith angles (z 0 = 85 • ).

This comparison allows us to say that this additional correction is not insignificant. Although to have an idea of the error on the estimated altitude of the object, it is sufficient to look at the Fig. 9 of the evolution of the lateral shift as a function of the zenith angle. In consistency with our intuition, the error increases with the zenith angle, it is about 10 m at z 0 = 60 • and a little less than 100 m at z 0 = 80 • . difference tells us about the additional directions available when a phase variance budget has been previously set for the adaptive optics system.

From Fig. 16, we observe that the use of the numerical solution of the lateral shift does not induce much gain in sky coverage. In fact, the difference between the two sequences of zenith angles does not exceed 0.1 • up to 50 • of zenith angle.

Based on the informations provided by both Fig. 15 and Fig. 16, we can reasonably say that in the case of adaptive optics, the first order approximation is quite sufficient because there is little difference between the phase variance calculated with either one, and this is directly reflected i n t he s ky c overage gain which is small too.

Effect of the lateral shift on the observation of nearby objects

Sometimes one may want to observe objects less than 100 km away from Earth, it is the case in meteors observations which are very much popular at the moment because of the creation of national and international networks of cameras dedicated to photographing meteorite falls such as FRIPON (Fireball Recovery and Inter Planetary Observation Network) [START_REF] Colas | Proceedings of the International Meteor Conference[END_REF](Colas et al. , 2020)). Usually, the goal of those projects is to estimate the land-However, these large distances must be handled with care, on the one hand the error on the falling point is difficult to estimate because there are many parameters to take into account such as the wind, the falling direction, the observation points, etc. And on the other hand, current photography networks are unable to take this correction into account, due to the low resolution of the cameras used (less than 10 arc minutes).

This correction is thus to be included in the uncertainty on the vertical position of the object, and can be used later to improve and refine future fireball photography systems.

Conclusion

We first established the equations of refraction based on the tracing of the light path of a ray through the atmosphere. In order to solve these equations, we used a two-layer atmosphere model with a fixed temperature gradient and variable temperature, pressure and relative humidity at the observer level. This led to a first estimate of the lateral shift b 0 using a fourth order Runge-Kutta integration. In order to simplify the calculations and drop the integration along the light path, we derived then a secondorder approximation of the lateral shift. Furthermore, we have rewritten this approximation for the model of dry atmosphere, and this results in a compact form that depends only on the apparent zenith angle, on the pressure and temperature values at the observer level, and on the main moments of the atmosphere. Starting from the second order approximation, we have extracted two other estimators of the lateral shift, one of which was found to be the lateral shift given by [START_REF] Wallner | Imaging Through the Atmosphere[END_REF]. Table 1 summarizes the four ways presented in this paper to compute the lateral shift. The approximation b (3/2) offers the great advantage of being separable into two function of the wavelength lambda and the zenith angle z 0 .

From the integration of the equations, we also investigated the impact of pressure, temperature and observation wavelength. The main meteorological factor was found to be pressure while the variation as a function of wavelength, crucial for several applications, does not exceed 3% in the visible and near infrared spectrum.

In addition to this, we attempted to evaluate the relevance of the numerical computation of the lateral shift for adaptive optics. Despite the fact that Wallner's formula overestimates the values of the lateral shift, we realized that the numerical integration was of little interest for adaptive optics, as the first order approximation is sufficient to study the effects of chromatic shear.

Finally, using numerical integration of refraction equations, we evaluated the importance of the lateral shift in compensating the impact of refraction on the observation of nearby objects such as meteors. We have concluded that this correction is far from being insignificant, especially when the object is close to the ground where it reaches 30% of the angular refraction (at infinity).

  error . . . . . . . . . . . . . . . 7 3.2 Influence of temperature and pressure . . . . . . 7 3.3 Dependence on observation wavelength . . . . . 8 Application to meteor trajectography and adaptive optics 9 4.1 Relevance of the second order approximation of the lateral shift for adaptive optics . . . . . . . . 9 4.2 Effect of the lateral shift on the observation of nearby objects . . . . . . . . . . . . . . . . . . . 10 Conclusion 11

Fig. 1 .Fig. 2 .

 12 Fig. 1. Bending of light rays due to atmospheric refraction.

Fig. 4 .

 4 Fig. 4. Evolution of the three characteristic lengths as a function of the initial altitude.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Evolution of each of the three terms according to the zenith angle.

Fig. 7 .

 7 Fig. 7. Absolute error of the three estimators.

Fig. 8 .

 8 Fig. 8. Relative error of the three estimators.

Fig. 12 .Fig. 13 .

 1213 Fig. 12. Absolute deviation of the lateral shift from its value at λ = 550 nm as a function of wavelength, computed using Eqs. 16.

Fig. 15 .

 15 Fig.15. The ratio ∆b 0 /∆b(1) 0 as a function of zenith angle.

Fig. 17 .

 17 Fig. 17. Correction ratio σ/R a versus object's altitude (H).

Fig. 16 .

 16 Fig. 16. The improvement in sky coverage induced by the use of the numerical solution instead of the first-order approximation, with λ SCI = 1630 nm.

Table 1 .

 1 Overview of the four estimators of the lateral shift.

	Estimator	Equation	Zenith angles where (X) r ≤ 1%.
	b 0 (integration of the
		system of Eqs. 16)
	b (1)	1	0 -55 •
	b (3/2)	48	0 -70 •
	b (2)	39	0 -75 •

Appendix A: A model for the dry atmosphere Below we present the derivation that leads to the expression of the optical index and its gradient, in the case of dry air and using the temperature profile presented in section 2.2. We begin by rewrinting the Gladstone-Dale relation:

Under the ideal gas approximation, and by decoupling the dependencies in ρ and λ (respectively density and wavelength), Eq. A.1 reads:

where A D is called the reduced refractivity of dry air and is given by several models, one of the most recent expressions being given by the equation 2 of Ciddor (1996). That is:

This expression is given for a wavelength in µm and A D is in hPa -1 K. The previous expression of A D applies for wavelengths from 300 nm to 1690 nm and is also valid over the near infrared spectrum. We then write P/T as a function of altitude h assuming a constant and known temperature gradient in the troposphere ω and nil temperature gradient above the tropopause. Hence, we have:

H t is the altitude of the tropopause. In this framework, the pressure at a given altitude h ≤ H t depends on the temperature as follows: -M D is the molar mass of dry air, -ω := dT /dh = -6.5 K/km is the temperature gradient in the troposphere, -R is the universal gas constant.

While in the stratosphere, assuming a constant temperature, one gets:

where T t is the temperature at the tropopause. We take H t = 11 km as in [START_REF] Hohenkerk | The Computation of an Angular Atmospheric Refraction at Large Zenith Angles[END_REF].

Hence, the refractive index with respect to height h reads : And its derivative is given by:

(A.9) Figure A.1 shows the evolution of the refractive index and its derivative as a function of altitude. We notice that n decreases exponentially in the stratosphere. Also, we observe a discontinuity of dn/dh at the boundary between the troposphere and the stratosphere, which is due to the discontinuity of the temperature gradient.


Jaime Arias
email: arias@lipn.univ-paris13.fr

Wojciech Penczek
email: penczek@ipipan.waw.pl

Laure Petrucci
email: petrucci@lipn.univ-paris13.fr

Teofil Sidoruk
email: t.sidoruk@ipipan.waw.pl

ADT2AMAS: Managing Agents in Attack-Defence Scenarios Demonstration Track

Keywords: asynchronous multi-agent systems, scheduling, attack-defence trees

Expressing attack-defence trees (ADTrees) in a multi-agent setting allows for studying a new aspect of security scenarios, namely how the number of agents and their task assignment impact the performance of attacking and defending strategies executed by agent coalitions. Our tool ADT2AMAS allows for transforming ADTrees into extended asynchronous multi-agent systems and computing an optimal schedule with the minimal number of agents. ADT2AMAS is integrated within the graphical verification platform CosyVerif, but can also be run standalone.

INTRODUCTION AND MOTIVATIONS

Security of safety-critical multi-agent systems [START_REF] Michael | An Introduction to Multiagent Systems[END_REF] is one of the major challenges. Attack-defence trees (ADTrees) [START_REF] Aslanyan | Pareto Efficient Solutions of Attack-Defence Trees[END_REF][START_REF] Kordy | Foundations of Attack-Defense Trees[END_REF] have been developed to evaluate the safety of systems and to study interactions between attacker and defender parties. They provide a simple graphical formalism, where nodes represent possible attacks against a system and defences employed to protect it. However, the defenders might forgo a costly defence, or a defence could take too much time and thus fail. The attack to be achieved is the root one, which requires all (for the AND nodes) or part (for the OR nodes) of those in its subtrees to be successful.

Recently, it has been proposed to model ADTrees in the formalism of asynchronous multi-agent systems (AMAS) extended with certain ADTree characteristics [START_REF] Arias | Hackers vs. Security: Attack-Defence Trees as Asynchronous Multi-agent Systems[END_REF][START_REF] Petrucci | Squeezing State Spaces of (Attack-Defence) Trees[END_REF]. In this setting, one can reason about attack/defence scenarios considering agent distributions over the tree nodes and their impact on the feasibility and performance (quantified by metrics such as time and cost) of attacking and defending strategies executed by specific agent coalitions. Both the number of agents available and their distribution over the ADTree nodes affect these performance metrics. Hence, there arises the problem of an optimal scheduling, that is, obtaining an assignment that guarantees the lowest possible execution time while using the minimum number of agents required for an attack to be feasible. The tool ADT2AMAS implements the algorithm of [START_REF] Arias | Minimal Schedule with Minimal Number of Agents in Attack-Defence Trees[END_REF] that finds such an assignment for a given ADTree. As the ADTree formalism introduces unique caveats, our approach differs from that of classical process scheduling [START_REF] Kwok | Static Scheduling Algorithms for Allocating Directed Task Graphs to Multiprocessors[END_REF][START_REF] Lombardi | Stochastic allocation and scheduling for conditional task graphs in multi-processor systems-on-chip[END_REF][START_REF] Singh | A Survey of Static Scheduling Algorithm for Distributed Computing System[END_REF]; to the best of our knowledge, this is the first work dealing with agents in this context. Notably, the This work is supported by the IEA project PARTIES and FNR/NCBiR project STV.

scheduling algorithm optimises both the number of agents and the attack time, and runs in quadratic time wrt. the number of nodes. However, it is applied to a number of models that is exponential in the number of OR and defence nodes.

FORMAL BACKGROUND

A well-known formalism [START_REF] Buldas | Rational Choice of Security Measures Via Multi-parameter Attack Trees[END_REF][START_REF] Fila | Exploiting Attack-Defense Trees to Find an Optimal Set of Countermeasures[END_REF][START_REF] Jhawar | Attack Trees with Sequential Conjunction[END_REF][START_REF] Kordy | DAGbased Attack and Defense Modeling: Don't Miss the Forest for the Attack Trees[END_REF][START_REF] Kumar | Effective Analysis of Attack Trees: A Model-Driven Approach[END_REF][START_REF] Mauw | Foundations of Attack Trees[END_REF][START_REF] Slater | Toward a Secure System Engineering Methodolgy[END_REF], ADTrees model security scenarios as an interplay between two opposing parties. Figure 1 depicts basic constructs, [START_REF] Arias | Hackers vs. Security: Attack-Defence Trees as Asynchronous Multi-agent Systems[END_REF] gives a comprehensive overview.

Attacking and defending actions are depicted in red and green, respectively. The leaves represent individual actions at the highest level of granularity. Different types of gates (available in both attack or defence types) allow for modeling increasingly broad intermediary goals, all the way up to the root, which corresponds to the overall objective. OR and AND gates are defined analogously to their logical counterparts. SAND is a sequential variant of the latter, i.e. the entire subtree 𝑎 𝑖 needs to be completed before handling 𝑎 𝑖+1 . Countering actions can be expressed using gates CAND (counter defence; successful iff 𝑎 succeeds and 𝑑 fails), NODEF (no defence; successful iff either 𝑎 succeeds or 𝑑 fails), and SCAND (failed reactive defence; sequential variant of CAND, where 𝑎 occurs first). ADTree nodes may have numerical attributes, e.g. the time needed to perform an attack, or its financial cost. Conditions are boolean functions over the attributes, used as constraints to counter-defence nodes. The translation of ADTrees to extended asynchronous multiagent systems (EAMAS, [START_REF] Arias | Hackers vs. Security: Attack-Defence Trees as Asynchronous Multi-agent Systems[END_REF][START_REF] Jamroga | Towards Partial Order Reductions for Strategic Ability[END_REF][START_REF] Jamroga | Towards Partial Order Reductions for Strategic Ability[END_REF]), called also models, was proposed in [START_REF] Arias | Hackers vs. Security: Attack-Defence Trees as Asynchronous Multi-agent Systems[END_REF][START_REF] Petrucci | Squeezing State Spaces of (Attack-Defence) Trees[END_REF]. To that end, the semantics of an asynchronous MAS of [START_REF] Jamroga | Towards Partial Order Reductions for Strategic Ability[END_REF][START_REF] Jamroga | Towards Partial Order Reductions for Strategic Ability[END_REF] was extended to account for the node attributes and conditional constraints. Each ADTree node corresponds to a separate automaton, with specific patterns used for different types of nodes. Transitions are decorated with actions and the following elements: (i) a message 𝑓 𝑀 (𝑡) ∈ ({!, ?} × 𝑀) ∪ {⊥} indicating whether there is a synchronisation on a message or the action is local; (ii) a guard In this setting, groups of agents working for the attacking and defending parties can be considered. Note that the feasibility of an attack is not affected by the number or distribution of agents over ADTree nodes, while some performance metrics, such as time, clearly depend on both factors (e.g. a lone agent can sequentially handle all the actions, albeit usually much slower). Furthermore, the success or failure of particular defence nodes results in a number of unique configurations, leading to multiple model variants to be analysed. Consequently, the optimal distribution of agent coalitions is of vital importance for both parties, allowing them to prepare for multiple scenarios, depending on how many agents they can afford to recruit. Thus, synthesising an assignment that achieves a minimal execution time using the least possible number of agents is a problem of high interest and importance.

𝑎 (a) leaf (attack) 𝐴 𝑎1 𝑎𝑛 • • • (b) AND 𝐴 𝑎1 𝑎𝑛 • • • (c) OR 𝐴 𝑎1 𝑎𝑛 • • • (d) SAND 𝑑 (e) leaf (defence) 𝐴 𝑎 𝑑 (f) CAND 𝐴 𝑎 𝑑 (g) NODEF 𝐴 𝑎 𝑑 (h) SCAND

ARCHITECTURE AND TECHNOLOGY

ADT2AMAS [START_REF] Arias | ADT2AMAS[END_REF] is an open source tool written in C++17. It allows for:

(1) translating an ADTree into an EAMAS, and (2) computing a minimal scheduling with a minimal number of agents to make the root node's task successful. The ADTree can be specified using the provided API, or loaded from a simple-syntax text file (Figure 2a). For (1), the tool generates a L A T E X file of the EAMAS, and an .imi file containing its specification, which can be read by the modelchecker IMITATOR [START_REF]IMITATOR[END_REF] in order to run automatic verifications. For (2), ADT2AMAS generates a L A T E X file of each tree transformation needed to compute the minimal assignment, and an ASCII table of the minimal scheduling. Compiling L A T E X to PDF provides a visual feedback of the EAMAS and of the scheduling algorithm's steps.

ADT2AMAS is also integrated within the CosyVerif formal specification and verification environment [START_REF]CosyVerif[END_REF], thus providing a multiplatform, user-friendly, zero-configuration tool (Figure 2b). Its Java-Script web interface, CosyDraw, allows for designing formal models extensible via FML (Formalism Markup Language) and GrML (Graph Markup Language) specifications. Furthermore, it supports any tool providing a SOAP web service. In our case, we extended the graphical interface with the ADTree elements (Figure 2c) and implemented a web service offering the algorithms supported by ADT2AMAS [START_REF] Arias | ADT2AMAS Alligator Service[END_REF]. The tool can be accessed at https://cosyverif.lipn.univ-paris13.fr, and its video demonstration at https://youtu.be/DGLtUSP-ao8.

EXPERIMENTS

Selected experimental results are provided in Table 1, with additional ones available on the tool web page [START_REF]ADT2AMAS Minimal Scheduling[END_REF]. Subsequent columns denote: (1) the case study, (2) the number of possible defence configurations, which influence the number of models to consider, (3) the number of a model variant; note that some combinations of defences result in the same model, or do not allow for a successful attack, (4) the minimal number of agents obtained by our scheduling algorithm, (5) the minimal attack time. The tool also provides a corresponding schedule of agents' actions (not shown here), and generates the resulting EAMAS (currently, for 1 agent only), which can then be passed to the IMITATOR model checker. This allows for the verification of other desired properties (e.g. a model with a higher attack time might be preferable due to its far lower cost).

FUTURE WORK

Currently, model specifications (.imi files) including more than one agent are prepared by hand. We plan to improve the integration of ADT2AMAS with IMITATOR by a fully automatic generation of EAMAS for an arbitrary number of agents. Furthermore, ADT2AMAS can be integrated with other model checkers such as UPPAAL [START_REF] Behrmann | UPPAAL 4.0[END_REF].

Figure 1 :

 1 Figure 1: Basic ADTree gates (or nodes)

Figure 2 :

 2 Figure 2: Graphical Interface for ADT2AMAS𝑓 𝐺 (𝑡) ∈ Guards constraining the transition; (iii) an update function 𝑓 𝑡 : AT 𝑡 → EXP (AT, FUN ) expressing how the transition modifies the attributes. In this setting, groups of agents working for the attacking and defending parties can be considered. Note that the feasibility of an attack is not affected by the number or distribution of agents over ADTree nodes, while some performance metrics, such as time, clearly depend on both factors (e.g. a lone agent can sequentially handle all the actions, albeit usually much slower). Furthermore, the success or failure of particular defence nodes results in a number of unique configurations, leading to multiple model variants to be analysed. Consequently, the optimal distribution of agent coalitions is of vital importance for both parties, allowing them to prepare for multiple scenarios, depending on how many agents they can afford to recruit. Thus, synthesising an assignment that achieves a minimal execution time using the least possible number of agents is a problem of high interest and importance.

Table 1 :

 1 Experimental results on case studies case study # def. no. of model # agents time

	treasure	2	1	2	125 min
	forestall		1	1	43 days
		4	2	1	54 days
			3	1	53 days
	iot-dev	4	1	2	694 min
	gain-admin		1	1	2942 min
		16	2	1	4320 min
			3	1	5762 min
	scaling	1	1	6	5 min


Julie Ledoux
email: julie.ledoux@ens-paris-saclay.fr

Luba Tchertanov

Receptor Tyrosine Kinase KIT: A New Look for an Old Receptor

Keywords: Receptor Tyrosine Kinase, RTK, full-length KIT cytoplasmic region, intrinsically disordered regions, phosphotyrosine, modelling, molecular dynamics, conformational transition, allosteric regulation, free energy landscape

RTK KIT regulates a variety of crucial cellular processes via its cytoplasmic (CD) domain composed of the tyrosine kinase domain crowned by highly flexible domains -juxtamembrane region, kinase insert domain and C-tail, key recruitment regions for downstream signalling proteins. We generated the 3D model of the full-length CD attached to the transmembrane helix to prepare a structural basis for characterization of the interactions of KIT with its signalling proteins (KIT INTERACTOME). This generic model of KIT in inactive state was studied by molecular dynamics simulation under conditions mimicking the natural environment of KIT. With the accurate atomistic description of the multidomain KIT dynamics, we explained its intrinsic (intradomain) and extrinsic (inter-domain) disorder, and represented the conformational ensemble of KIT through free energy landscapes. Strongly coupled movements within each domain and between distant domains of KIT prove the functional interdependence of these regions, described as allosteric regulation, a phenomenon widely observed in many proteins. We suggested that KIT in inactive state encodes all properties of the active protein and post-transduction events.

Receptor tyrosine kinases (RTKs) control various signalling pathways in cells. Their remarkable conformational plasticity enables the specific recognition of many molecules such as ligands, substrates or proteins. In solution, RTKs are at equilibrium between different conformations ranging from an inactive autoinhibited state to a fully active state. Ligand-induced activation of RTKs leads to recruitment and activation of multiple downstream signalling proteins which alter the expression of genes governing cell physiology [START_REF] Lemmon | Cell Signaling by Receptor Tyrosine Kinases[END_REF]. Explicit elucidation of signalling events is an important and unresolved problem in cell biology. The initiation of these cascade-like processes involves different domains of RTK, each performs specific actions, finely concerted by a regulated allosteric mechanism controlling all functional biological processes [START_REF] Gunasekaran | Is allostery an intrinsic property of all dynamic proteins?[END_REF][START_REF] Wodak | Allostery in Its Many Disguises: From Theory to Applications[END_REF].

Focusing on the RTK KIT, an important target in oncology [START_REF] Pham | KIT and Melanoma: Biological Insights and Clinical Implications[END_REF][START_REF] Gilreath | Novel approaches to treating advanced systemic mastocytosis[END_REF], we will discuss this RTK as a key regulator of intracellular signalling mediated by regions possessing multiple phosphorylation sites: juxtamembrane region (JMR), kinase insert domain (KID), activation loop (A-loop) and C-terminal tail. Since these regions are very flexible or disordered, their properties are not yet well understood. As KIT is a multi-functional protein, its different regions regulate catalytic processes and/or events that activate and control the signalling cascade. To complete the functions required in more than one region, these regions should be directly or collaterally coupled.

We generated the 3D model of the inactive full-length KIT attached to the transmembrane helix (Fig. 1.), to (i) prepare a structural basis for the characterisation of interconnections between functional regions -JMR, tyrosine kinase (TK) domain, KID and C-tail -, and (ii) establish the interactions of KIT with its signalling proteins. Then, we investigated this model by molecular dynamics (MD) simulation in conditions mimicking its natural environment. We suggested that such atomistic description of KIT will fully elucidate structural and dynamical properties of its different functional regions. To the best of our knowledge, we have presented for the first time a model of a full-length cytoplasmic region of an RTK KIT attached to a transmembrane helix and its molecular dynamics simulations under conditions that mimic its natural environment [START_REF] Ledoux | The Inherent Coupling of Intrinsically Disordered Regions in the Multidomain Receptor Tyrosine Kinase KIT[END_REF]. Analysis of the simulation data (three 2-µs MD trajectories) put in evidence that the multidomain RTK KIT is a modular protein consisting of a quasistable TK domain crowned by at least four intrinsically disordered (ID) regions -JMR, KID, A-loop and C-tail. These ID regions belong to two types -the very elongated (extended) and poorly folded regions (JMR, A-loop and C-tail), and the globule-like (collapsed) KID having a high level of the helical structures. KIT ID regions contain transient structures (helices and β-strands) and their local architecture displays various degrees of compaction and elongation. Therefore, the structure of each ID domain of KIT represents a very complex mixture of a broad variety of differently folded conformations which describe a reversible folding-unfolding process, specific for each ID domain. In particular, the KID, composed of transient helices linked by coils, is the most disordered domain in respect to other KIT domains, but shows a globule-like shape stabilised by noncovalent interactions [START_REF] Inizan | The First 3D Model of the Full-Length KIT Cytoplasmic Domain Reveals a New Look for an Old Receptor[END_REF][START_REF] Ledoux | Folding and Intrinsic Disorder of the Receptor Tyrosine Kinase KIT Insert Domain Seen by Conventional Molecular Dynamics Simulations[END_REF]. Also, the inherently disordered KID shows different positions derived from two types of motions -linear (translation) and angular (rotation) displacement -regarding the stable TK domain. The elongated regions (JMR, A-loop and C-tail) show rather local disorder as evidenced by alternating positions of their short segments relative to the stable TK domain.

The two-level disorder (intrinsic and extrinsic) provides high conformational variability of KIT and supplies the high adaptability of JMR, KID and C-tail required for the scaffolding (docking sites) and recruitment of different protein partners of KIT, and accomplishes the tight regulation of cellular processes. Consequently, the overall structure of KIT represents a continuous spectrum of conformations with a different degree and depth of disorder, thereby generating a complex protein structural space. It is partially reflected by free energy landscapes lacking a unique global deep minimum as typically observed in ordered proteins. Such energy landscapes, with two local minima joined by a 'flattened plateau' containing the intermediate conformations, show that KIT is extremely sensitive to different environmental changes (e.g. phosphorylation) that can alter its free energy landscape in different ways. We suggest that JMR, displaced from its packed auto-inhibited position upon the SCF-induced activation of KIT, will achieve higher levels of disorder, and therefore a higher level of adaptability for the recruitment of signalling proteins.

Since ID domains are multiple in RTK KIT, does the disorder/order of one domain depend on the disorder/order of other remote regions? As evidenced by the cross-correlations, the highly coupled motions of distant regions of KIT suggest their functional dependence which is classified as allosteric regulation, phenomenon largely observed in many proteins. In particular, the coupling motions within the TK domain reflect the allosteric regulation of kinase function which is well-described for different non-receptor and receptor tyrosine kinases [START_REF] Ahuja | Dynamic allostery-based molecular workings of kinase:peptide complexes[END_REF]. The coupled/uncoupled motions of A-loop and JMR were described through their allosteric communication in the wild-type KIT and oncogenic mutants [START_REF] Laine | Allosteric Communication across the Native and Mutated KIT Receptor Tyrosine Kinase[END_REF][START_REF] Allain | Allosteric pathway identification through network analysis: from molecular dynamics simulations to interactive 2D and 3D graphs[END_REF].

Fig. 1 .

 1 Fig. 1. The 3D model of the inactive full-length KIT attached to the transmembrane helix. Modelling and study by MD simulations (top panel) analysed to (i) characterise the global motions and coupling, (ii) to represent the MD conformations as the free energy landscape, (iii) to describe the relative positions of the KIT functional regions and the position of phosphotyrosine residues.