Jessica Colombel

David Daney

François Charpillet

Holistic view of Inverse Optimal Control by introducing projections on singularity curves

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Holistic vision of Inverse Optimal Control

Jessica Colombel, David Daney, Franc ¸ois Charpillet Abstract-Inverse Optimal Control (IOC) is a problem that exists in many fields. In the context of human motion analysis, many methods for solving this problem have been proposed. This paper presents the Projected Inverse Optimal Control (PIOC) approach which puts forward a simple and complete vision of the problem. PIOC is composed of several independent elements allowing to cover problems of dimensioning, measurement noise and reliability of solution. The first element is the parameterization of the data, allowing to reduce the dimensions of the problem. The second element corresponds to the decoupling of the conditions including the conditions of identifiability of the parameters, singularity of the problem and feasibility of the solution. The last element refers to the choice of solution. We show with this approach that the classical solution methods are in fact projections in the space of trajectories. PIOC also allows us to propose a simple algorithm for the choice of basis, a recurrent problem in IOC.

Inverse Optimal Control, Human Motion Analysis, Identification.

I. INTRODUCTION

The Inverse Optimal Control (IOC) problem is a complex problem studied in many disciplines. In particular, the IOC is related to the Inverse Reinforcement Learning problem whose expression and proposed resolution methods differ [START_REF] Azar | From inverse optimal control to inverse reinforcement learning: A historical review[END_REF]. For human motion analysis, the IOC is a well known problem. Thanks to this method, the literature analyzes the movement in order to imitate it with a robot [START_REF] Mombaur | From human to humanoid locomotion-an inverse optimal control approach[END_REF], [START_REF] Park | Inverse Optimal Control for Humanoid Locomotion[END_REF], to predict human movements to better interact with humans [START_REF] Mainprice | Predicting human reaching motion in collaborative tasks using Inverse Optimal Control and iterative re-planning[END_REF], [START_REF] Gaurav | Discriminatively Learning Inverse Optimal Control Models for Predicting Human Intentions[END_REF] or even to split the movement into subtasks [START_REF] -S. Lin | Human motion segmentation using cost weights recovered from inverse optimal control[END_REF]. Two methods have been particularly used in the literature: the Bilevel method [START_REF] Mombaur | From human to humanoid locomotion-an inverse optimal control approach[END_REF], [START_REF] Albrecht | Modeling and Analysis of Human Navigation with Crossing Interferer Using Inverse Optimal Control[END_REF] and the AIOC method [START_REF] Panchea | Inverse Optimal Control for Redundant Systems of Biological Motion[END_REF], [START_REF] Puydupin-Jamin | A convex approach to inverse optimal control and its application to modeling human locomotion[END_REF]. In the context of human motion analysis, we are particularly interested in IOC under uncertainties: noisy data, uncertain model, unknown cost function and parameter selection.

This paper proposes a new approach for Inverse Optimal Control problem, when observed data are noisy, which raise up the problem of optimality of the observed trajectory. We called it Projected Inverse Optimal Control (PIOC). This approach is divided into three elements, the parametrization of the measurement, the decoupling of conditions of viability of the problem and the choice of solution. Our approach brings a simple and complete vision of the IOC resolution problem. In particular, we show that the resolutions methods proposed in the literature are equivalent to particular projections in the trajectory space. The choice of the resolution method is equivalent to a choice of solution according to different criteria.

Thus, the originality of this paper lies in a new simple vision of the IOC problem through:

• parameterization of the measurement states and its impact on the dimensionality and identifiability of the problem;

• decoupling of the conditions: identifiability, singularity and feasibility conditions; • illustration of the choice of the solution by a projection;

• solution for the choice of the cost function. First, this paper introduces the problem statement of IOC by defining Direct Optimal Control (DOC) and Bilevel and AIOC resolution methods. Then, the PIOC approach is presented by detailing its different steps: parametrization, decoupling conditions and choice of solution. We also provide a simple example to illustrate our approach. Finally, we will discuss the perspectives and limitations of the proposed approach before concluding this paper.

II. PROBLEM STATEMENT

A. Context

The analysis of human movement by the inverse optimal control is based on the assumption that the movement is optimal according to a criterion. Bipedal walking is considered as an optimal movement [START_REF] Alexander | The Gaits of Bipedal and Quadrupedal Animals[END_REF], as well as simpler movements such as picking where the hand is following a trajectory minimizing the criterion of torque variation [START_REF] Uno | Formation and control of optimal trajectory in human multijoint arm movement[END_REF]. Moreover, the biological sensorimotor control would itself be governed by optimal laws [START_REF] Todorov | Optimality principles in sensorimotor control[END_REF]. Given the multiplicity of movements and their representation, many criteria have been proposed to generate or analyze through IOC human movement [START_REF] -S. Lin | Human motion segmentation using cost weights recovered from inverse optimal control[END_REF], [START_REF] Panchea | Inverse Optimal Control for Redundant Systems of Biological Motion[END_REF], [START_REF] Carreno-Medrano | Analysis of Affective Human Motion During Functional Task Performance: An Inverse Optimal Control Approach[END_REF]- [START_REF] Rebula | A Robustness Analysis of Inverse Optimal Control of Bipedal Walking[END_REF]. These criteria, or cost functions, are the basis for the expression of DOC. While the DOC allows to generate trajectories, the IOC allows to analyze trajectories and to understand how they were generated.

The following subsections will present the DOC as well as two methods for solving the IOC problem. We will simplify their presentation as much as possible in order to stay pedagogical. Finally, we will discuss the limitations of these methods.

B. Direct Optimal Control

In the context of human (or robot) motion, the aim of the DOC is to generate a trajectory s * optimal for a convex criterion

C(s) = nc k ω k .C k , composed by a sum of weighed (ω = [ω 1 , . . . , ω nc ]) cost functions (C k , k = 1, . . . , n c ), under n f equality constraints f i (s) (i = 1, . . . , n f ) and n h inequality constraints h j (s) (j = 1, . . . , n h ). s * = arg min s C(s) s.t. f 1,...,n f (s) = 0, h 1,...,n h (s) ≤ 0 (1)
where s ∈ R ns×nt is a trajectory composed, by n t of frames; each are defined by a state s t of dimension n s , such that s = [s 1 , . . . , s t , . . . , s nt ] T .

The existance of a solution to this problem is assured by the convexity of the base [START_REF] Boyd | Convex Optimization[END_REF]. For that, the litterature suggest convexed cost functions (e.g. quadratic) and positive weights ω [START_REF] Panchea | Inverse Optimal Control for Redundant Systems of Biological Motion[END_REF].

C. Inverse Optimal Control

The IOC is the problem of finding the unknown weights ω associated to known cost functions. The constraints are also known. A lot of resolution methods were proposed to find the solution of the problem of IOC in the context of human movement analysis ( [START_REF] Mombaur | From human to humanoid locomotion-an inverse optimal control approach[END_REF], [START_REF] Park | Inverse Optimal Control for Humanoid Locomotion[END_REF], [START_REF] Clever | An Inverse Optimal Control Approach for the Transfer of Human Walking Motions in Constrained Environment to Humanoid Robots[END_REF], [START_REF] Berret | Evidence for Composite Cost Functions in Arm Movement Planning: An Inverse Optimal Control Approach[END_REF], [START_REF] Jin | Inverse Optimal Control with Incomplete Observations[END_REF]) but we will focus on two of them: the Bilevel resolution method and the Approximate IOC method.

1) Bilevel IOC: In the context of human motion analysis, the Bilevel IOC resolution method has been used for locomotion [START_REF] Mombaur | From human to humanoid locomotion-an inverse optimal control approach[END_REF], [START_REF] Clever | Inverse optimal control based identification of optimality criteria in whole-body human walking on level ground[END_REF] and for the arm [START_REF] Albrecht | A bilevel optimization approach to obtain optimal cost functions for human arm movements[END_REF]. It is a resolution method that uses two levels of calculation (lower-level: DOC that generates optimal trajectory; upper-level: distance of generated trajectory to observed trajectory) that loop until they converge to the result: the optimal trajectory closest to the observed one [START_REF] Sinha | A Review on Bilevel Optimization: From Classical to Evolutionary Approaches and Applications[END_REF]. The IOC Bilevel is expressed as follows:

min ω ∥s * -s M ∥ 2 with s * = argmin s nc k=1 ω k nsnt t=1 C k (s t ) s.t. f 1,...,n f (s) = 0, h 1,...,n h (s) ≤ 0 (2)
with s M the observed trajectory. We can already note in this expression that if the data are noisy, we converge to the optimal trajectory closest to the noisy one.

2) Approximated IOC: The so-called "Approximate" resolution methods are methods based on properties of an optimal problem and resulting conditions. The best known conditions are the Karush-Kuhn-Tucker conditions. They have been used a lot in the framework of IOC [START_REF] -S. Lin | Human motion segmentation using cost weights recovered from inverse optimal control[END_REF], [START_REF] Puydupin-Jamin | A convex approach to inverse optimal control and its application to modeling human locomotion[END_REF], [START_REF] Carreno-Medrano | Analysis of Affective Human Motion During Functional Task Performance: An Inverse Optimal Control Approach[END_REF], [START_REF] Panchea | Towards solving inverse optimal control in a bounded-error framework[END_REF]- [START_REF] Aghasadeghi | Inverse optimal control for differentially flat systems with application to locomotion modeling[END_REF]. Other conditions have also been used for IOC, such as the Jacobi-Bellman equations [START_REF] Pauwels | Inverse optimal control with polynomial optimization[END_REF]- [START_REF] Johnson | Inverse optimal control for deterministic continuous-time nonlinear systems[END_REF].

The KKT conditions are expressed as follows:

nc k=1 ω k ∂C k ∂s (s * ) + n f i=1 λ i ∂f i ∂s (s * ) + n h j=1 ν j ∂h j ∂s (s * ) = 0 stationnary f i (s * ) = 0, i = 1, . . . , n f h j (s * ) ≤ 0, j = 1, . . . , n h primal feasibility
ν j ≥ 0, j = 1, . . . , n h dual feasibility ν j h j (s * ) = 0, j = 1, . . . , n h complementary slackness

(3) The stationnary condition can also be written with matrices:

[J ω , J λ , J ν ] J   ω λ ν   z = 0 (4)
with J the Identification Matrix and J ω , J λ and J ν the sub-matrices associated with the parameters ω (what the IOC wants to retrieve) and λ, ν the Lagrangian coefficients of the constraints of equality and inequality, respectively.

The basis C is convex, so that these conditions are considered necessary but also sufficient for the optimality of the problem [START_REF] Hanson | On sufficiency of the Kuhn-Tucker conditions[END_REF].

In the literature, these conditions are used to solve the IOC problem by using them as "residuals". The solution is then done by minimizing residuals [START_REF] -S. Lin | Human motion segmentation using cost weights recovered from inverse optimal control[END_REF], [START_REF] Puydupin-Jamin | A convex approach to inverse optimal control and its application to modeling human locomotion[END_REF], [START_REF] Carreno-Medrano | Analysis of Affective Human Motion During Functional Task Performance: An Inverse Optimal Control Approach[END_REF], [START_REF] Panchea | Towards solving inverse optimal control in a bounded-error framework[END_REF].

D. Limitations of the methods

Each of these resolution methods has limitations that do not allow for optimal use in human motion analysis. The main drawback of Bilevel is the computation time for the double optimization loop [START_REF] Johnson | Inverse optimal control for deterministic continuous-time nonlinear systems[END_REF]. The approximate methods pose problems of robustness [START_REF] Aswani | Inverse Optimization with Noisy Data[END_REF]. In particular, it is necessary to know to what extent one can trust the results obtained [START_REF] Colombel | On the Reliability of Inverse Optimal Control[END_REF].

However classical IOC methods work well if the problem is well posed and the data are perfect. In particular, if the basis chosen for the IOC correspond well to the one used to generate the reference trajectory, which must be provided without uncertainties, the resolution methods previously mentioned will find the solution without difficulty. In practice, this case cannot be considered: it is then necessary to take into account the natural uncertainty of the problem, in particular for the analysis of human motion which interests us particularly.

Moreover according to Ab Azar [START_REF] Azar | From inverse optimal control to inverse reinforcement learning: A historical review[END_REF], the IOC problem suffers more broadly from several problems:

• Ill-posedness problem (existence and uniqueness of the solution); • Non-convexity of functions and/or constraints;

• Data availability (incomplete or noisy data); • Non-linearity, complexity and dimensionality;

• Feature selection (relevance of parameters for dimensionality reduction); • Generalization (choice of cost functions/results obtained). In our approach, we want to address these difficulties and propose elements to overcome them. We therefore propose to unify the different methods by highlighting the central role of the projection of uncertain trajectory on an ideal provided by the choice of bases.

III. PSEUDO INVERSE OPTIMAL CONTROL APPROACH

Our approach, called PIOC, is is particularly interesting in the presence of uncertainties, whether they are on the data or the choice of cost functions or constraints.

A. Illustration

This article corresponds to the presentation of our general vision of the IOC. Thus we illustrate the different notions of our approach from a simple example already used in litterature [START_REF] Jin | Inverse Optimal Control with Incomplete Observations[END_REF]. We propose the same example in our previous article were the detailed equations are given [START_REF] Colombel | On the Reliability of Inverse Optimal Control[END_REF]: a 2 bars with two degrees of freedom θ 1 et θ 2 . Both the DOC and the IOC will be modeled from n c = 2 cost functions of torques:

C = τ 2 1 , τ 2 2 .
For each degree of freedom, the state is represented by the angle, its velocity and its acceleration. There is no inequality constraint in this example.

B. Parametrization: reduction of the dimension

Parameterization consists in expressing the observed state by parameters such as polynomials or Fourier series. Among the IOC problems cited earlier II-D, trajectory parameterization helps with some of them.

The parameterization allows to reduce the observed parameters, leading to a reduction of the dimensions of the problem while keeping all the information describing the state of the system. This reduction is necessary because the dimension of the state variables can be very large for human motion. For instance, if we observed a human arm we can have 8 angles and their velocity and acceleration, wich lead to n s = 24 state variables. Now suppose there are 30s of measurement at 25Hz, so n t = 750, the number of state variables will be n s * n t = 18000. By reducing the parameters of the problem, the time to solve the problem is reduced accordingly.

The data availability problem is also addressed by the parametrization. The data needs to be interpolated in order to be parameterized. This not only reduces the effects of missing data, but also reduces the impact of noise for the derived data thanks to the analytical expression of the parametrization.

For our approach, we propose to use a polynomial parameterization. This parametrization is already used in the literature to represent human motion trajectories [START_REF] Taylor | Polynomial approximations of gait for human motion analysis and visualization[END_REF], [START_REF] Hogan | Sensitivity of Smoothness Measures to Movement Duration, Amplitude, and Arrests[END_REF], and even for IOC [START_REF] -S. Lin | Human motion segmentation using cost weights recovered from inverse optimal control[END_REF], [START_REF] Panchea | Towards solving inverse optimal control in a bounded-error framework[END_REF]. Note that in our approach, we use the parametrization of the system states and not the cost functions, which is then called polynomial optimization [START_REF] Lasserre | Inverse polynomial optimization[END_REF].

For the observation of human motion, expressing the state by polynomials not only reduces the number of observed parameters, but also the number of constraints. The constraints concerning the start and end points and the kinematics are removed because they are directly included in the coefficients of the polynomial.

Let Q(t) be a polynomial describing the trajectory of a state parameter q, such that Q(t) = n d k=0 α k t k , with n d the degree of the polynomial. To simplify the calculations, we choose to define the starting and ending point at t = 0 and t = 1, respectively. It is then possible to associated start and end points with the polynomial coefficients and its derivation. The representation of the trajectories and their derivatives (first and second) by polynomials can be represented by a linear system in polynomial coefficients. The number of rows and columns of this linear system correspond to the initial and final constraints and to the number of parameters parameterizing the trajectory, respectively. The difference between the two gives the dimension of the space defining the set of trajectory allowing to go from an initial point to a final point. Theoretically, the larger this space is, the finer it can represent the observed trajectories.

If n d = 6, we obtain:

       
1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 2 0 0 0 0 1 1 1 1 1 1 1 0 1 2 3 4 5 6 0 0 2 6 12 20 30

                  α 0 α 1 α 2 α 3 α 4 α 5 α 6           -         q 0 q0 q0 q 1 q1 q1         = 0 (5)
This system is made of 6 equations for 7 unknowns, so it is possible to define 6 of the parameters in function of the 7th one. We can then express the polynomial Q only in terms of α 6 , as follows: (6) with a, b, c, d, e and f scalars only depending on the starting and ending points.

Q(t) = a + b t + c t 2 + (d -α 6 ) t 3 + (e + 3 α 6 ) t 4 + (f -3 α 6 ) t 5 + α 6 t 6
In this example, α 6 is the free coefficient which can express the whole polynomial. Note that if n d = 7, there will be 2 free parameters for each state. We propose to use the free coefficients as parameters for the IOC problem. It is important to understand that if n d ≤ 5, there is no free parameter, and therefore there is only one trajectory that corresponds to the start and end points. If there is only one trajectory described by this constraints, this trajectory is not the result of an optimization process but simply of constraints. Therefore, it is not possible to use a polynomial degree strictly lower than 6 to describe a trajectory supposedly resulting from an optimization. We insist on this lower bound because parameterizations of human motion tend to 5 th -order polynomial [START_REF] Hogan | Sensitivity of Smoothness Measures to Movement Duration, Amplitude, and Arrests[END_REF] for the study of the Jerk, but in the context of this optimization problem, this degree is not sufficient.

C. Decoupling of conditions

In our approach, the main principle is the decoupling of the conditions. As already introduced in our previous paper [START_REF] Colombel | On the Reliability of Inverse Optimal Control[END_REF], the solution of the problem is based on the verification of different conditions derived from KKT conditions. We consider an acceptable result when each of these conditions is met:

• Identifiability conditions; • Singularity condition; • Feasibility conditions. 1) Identifiability conditions: For a problem such as IOC, to admit a solution, all parameters must be identifiable. As detailed in our previous paper [START_REF] Colombel | On the Reliability of Inverse Optimal Control[END_REF], sub matrices (J ω , J λ , J ν ) associated with parameters need to be full ranks. If not a QR decomposition is suggested to eliminate the superfluous parameters.

Dimensionality: The identifiability is also dependent on the dimensions of the problem. In particular, the problem must have more (or at least equal) parameters than equations (cost functions and constraints). Note that in the context of parametric approach, as we proposed, for the matrix, this translates into having as many or more rows (of size (n d -5) * n q ) than columns (of size n c + (n f -n f p ) + n h with n f p constraints eliminated by the parametric). In other words (n d -5) * n q ≥ n c + (n f -n f p ) + n h . This means that for a given problem, model and constraints impact directly the degree of the polynomial allowing to solve the problem.

For example with our illustration III-A: polynomial trajectories eliminate all the equality constraints (n f = n f p ) from our problem and there is no constraints of equality (n h = 0). If the polynomials are of degree n d = 6, there are only two free parameters (one for each angle) so the matrix obtained is of size 2 × 2. If we want to study n c = 4 cost functions, we must increase the degree of the polynomial to n d = 7. This link between the degree of the polynomial and the problem studied is very important to take into account. All the more so as one must be careful since the higher the degree of the polynomial, the more likely the trajectory will have unwanted curve oscillations. Fig. 1: Illustration of a singularity curve from the example in section III-A. Two axis represents the free parameters α 6 and β 6 from the angle θ 1 and θ 2 , respectively. The color map is the determinant of the matrix det(J T J).

2) Singularity condition: Based on the KKT stationnary condition, this condition is a necessary but not sufficiant condition to express the optimality of the trajectory. We call it singularity condition because the Identification Matrix J needs to be singular to admit a non-trivial solution.

The singularity of a matrix is equivalent to at least one of the singular values is null or the determinant of the matrix is null. For our problem, the matrix J is often non-square, so the determinant should be expressed as: det(J T J). The determinant has the advantage of having an analytical expression, unlike singular values.

With our example, it is possible to display this condition in the trajectory space. We call the curve obtained the "singularity curve" (Fig. 1).

Note that this determinant is a necessary but not sufficient condition (KKT conditions on constraints are necessary to be sufficient). Figure 1 proposes a zoom on a portion (red square) on which all the trajectories from the DOC are.

3) Feasibility conditions: In our approach, we distinguish two types of feasibility conditions: the conditions associated to the constraints and the condition associated to the convexity of the problem.

As the KKT conditions express it, the conditions associated with the constraints are feasibility conditions. In a context of uncertainty, the constraints may not be respected. We consider that the verification of these constraints should be done according to the desired accuracy. This verification can even be optional.

The condition associated to the convexity of the problem is defined by the DOC (sec. II-B): the positivity of the weights obtained are necessary. Thus if there are not positives, the results are not considered reliable. It is important to understand that in our approach it is a condition and not a constraint. For example, in the litterature, AIOC methods use the positivity of ω as a constraints [START_REF] -S. Lin | Human motion segmentation using cost weights recovered from inverse optimal control[END_REF], [START_REF] Panchea | Towards solving inverse optimal control in a bounded-error framework[END_REF], [START_REF] Englert | Inverse KKT: Learning cost functions of manipulation tasks from demonstrations[END_REF]. We consider that the constraint forces the ω to be positive and thus modify the result to fit the constraint, even if it means a huge deviation from the observed trajectory. Using it as a condition allows us to invalidate the result if the weights are negative.

D. Choice of solution: projection

As previously mentioned, our approach is reserved for the problem of IOC under uncertainty. The parameterization and the decoupling of the conditions allow to tackle the difficulties of selection of the parameters, dimensionality, convexity and expression of the inverse problem. Thus, if the data are perfect, these steps are sufficient to obtain a reliable result whatever the chosen solution method. If the data are perfect and each of the previous conditions are met, the observed trajectory will be on a singularity curve. However, if the data are noisy the observed trajectory will not be on a singularity curve. Therefore, the resolution methods correspond to a projection of the observed trajectory to a singular trajectory on the singularity curve. Different methods do not project in the same way.

1) Type of projection: In this paper, we have mentioned two methods used in the literature for the application of IOC to human motion analysis: Bilevel and AIOC (Sec. II-C). These methods of solving the IOC problem correspond to different projections, one orthogonal and the other one that we will call "iso-weighted". Figure 2 illustrates these two projections for two uncertain trajectories.

a) Iso-weighted projection: The iso-weighted projection is obtained with the resolution of an unconstrained AIOC on ω. It corresponds to this equation:

z = arg min z ∥Jz∥ 2 s.t. det(J(s) T J(s)) = 0 (7)
Fig. 2: Comparison of projections for two trajectories in green and yellow on a singularity curve graph (Fig. 1). For each trajectory, the projection in plain line is iso-weighting and the dotted lines are "orthogonal". Note that for the green trajectory has two dotted projections, one has weight constraints (ω ≥ 0) the other has no constraints.

In reality, we only use a Singular Value Decomposition and the solution to this minimization should be the singular vector associated to the smallest singular value (more details in [START_REF] Colombel | On the Reliability of Inverse Optimal Control[END_REF]). Figure 3 shows for each trajectory, which weight ω is associated. In our example, ω is a vector of size n c = 2, so it is possible to express this vector as an angle. The yellow/red areas are the areas in which the ω is positive. They are therefore the areas in which it is possible to obtain a solution to the IOC problem with this resolution method.

It is to be noted that the resolution of AIOC is not a real "projection" on the trajectory space. Indeed, as the goal of the IOC is only to recover the omega, there is no real need to find the optimal trajectory corresponding to the weights found. We present it as a projection to show which singular trajectory the uncertain observed trajectory corresponds to. It is interesting to understand that for all iso-weighted trajectories, the minimum cost corresponding to this weighting is on the singularity curve. If the objective of the resolution is to minimize the cost of the observed trajectory, this solution is to be preferred.

b) Orthogonal Projection: The resolution of the IOC problem with a Bilevel corresponds an orthogonal projection constrained by positive omega. The high level is a minimization of the distance between two trajectories and in the polynomial framework the distance between two polynomials is linearly dependent on the free parameters. Let ξ p be the path error between two polynomials Q(t) et Q ′ (t), such that:

ξ p = 1 0 (Q(t) -Q ′ (t)) dt (8)
Fig. 3: Illustration of the iso-weighted area. Two axis represents the free parameters α 6 and β 6 from the angle θ 1 and θ 2 , respectively. The black curve is the singularity curve. The colormap is associated to the angle (in degree) that represents the vector ω of size 2. In yellow/red, ω is positiv.

If the constraints of the model and if the start and finish are identical, according to the equation (Eq. 6), we obtain:

ξ p = 1 0 (-(α 6 -α ′ 6 ) t 3 + 3 (α 6 -α ′ 6 ) t 4 -3 (α 6 -α ′ 6 ) t 5 + (α 6 -α ′ 6 ) t 6 ) dt = ∆ α 1 0 (-t 3 + 3 t 4 -3 t 5 + t 6 ) dt = - ∆ α 140 (9)
with ∆ α = α 6 -α ′ 6 . The same reasoning can be applied to each observed state as well as with a higher polynomial degree. We therefore have an error that is linear between the error on the trajectory and the error on the parameters. Now, in the space of parametric trajectories, the minimum squared norm is the shortest distance, which corresponds to an orthogonal projection between the measured trajectory parameters and those of the optimal trajectories. Moreover, since the low level is a DOC, ω are necessarily positive. Choosing the solution corresponding to a Bilevel then consists in taking the optimal trajectory closest to the observed trajectory. In other words, an orthogonal projection constrained by positive omega.

Orthogonal projection without constraints: As we said before, it does not seem relevant to us to constrain the positivity of ω but rather to use it as a condition (Sec III-C). This choice of projection which would be purely orthogonal would correspond to a method also used in the literature: the Onelevel IOC [START_REF] Hatz | Estimating Parameters in Optimal Control Problems[END_REF]. This method is composed by the high level of the bilevel but under KKT condition constraints. Onelevel IOC has the advantage of being faster than the Bilevel which has 2 levels of optimization. Moreover, according to our approach, we even propose to rewrite the Onelevel method as follows:

s * = arg min s ∥s -s M ∥ 2 s.t. det(J(s) T J(s)) = 0 (10)
As already mentioned in the paper, the use of the determinant rather than all the KKT conditions brings the advantage of being simpler and especially of decoupling the feasibility problem from the resolution, which is at the heart of our approach.

2) Choice of basis: The approach proposed in this paper aims to address many problems of IOC resolution in the context of human movement. The choice of the cost functions used to generate the observed trajectory is one of them. Figure 4 shows three singularity curves corresponding respectively to torque, angular power and acceleration multiplied by torque. If we observe a trajectory in the middle of its curves, how can we choose the best basis?

We propose a simple algorithm to choose a cost functions basis. This algorithm relies on properties of singularity curves and conditions proposed by our approach.

In a first step, the algorithm proposes to compute the determinants of each of the identification matrices associated to their basis. Once done, basis are sorted in ascending order for the determinants, so that the smallest determinant is first tested with a IOC resolution method. Doing the resolution of IOC on the smallest determinant allows the algorithm to be faster, especially if many basis are tested. However, this criterion cannot be unique because the ω also give clues on the reliability of the result. Thus, if ω is not positive, the next basis is tested. The first basis tested with a positive ω vector corresponds to the chosen basis. If no basis tested gives positive ω, then the algorithm returns that there are no basis associated with this trajectory.

Note that this algorithm is generic and can be used with any solving method (for example, AIOC or Bilevel). If Bilevel is chosen and the computation time of the algorithm is not a problem, it is possible to choose as selection criterion the basis corresponding to the optimal trajectory closest to the observed trajectory.

Be careful, the use of the determinant as a comparison criterion between different basis raises the question of the normalization of matrices.

IV. DISCUSSION AND PERSPECTIVES

The PIOC approach provides a holistic and generic vision of the IOC problem. In particular, it provides answers to issues of the IOC problem in the context of human motion: model, dimensionality, measurement noise, etc. Each of the proposed solutions is independent of the others and our approach is meant to be generalizable. We show in this paper that two methods of resolution of the IOC problem, Bilevel and AIOC, were just two different types of projection in the trajectory reference frame. This approach opens up many elements. Each of the steps proposed in our approach can be transposed with other methods. Concerning the parametrization, it would be possible to test with other parameters than the polynomial. Similarly, the optimality conditions on which we rely are the KKT conditions. The same approach can be taken with the Jacobi-Bellman conditions. As we work in an uncertain framework, the main interest of this approach is the decoupling between the resolution (or the choice of solution) and the different conditions allowing to choose if the result obtained is valid or not. This approach opens the door to many perspectives: better choice of solution according to the application, better choice of cost functions basis, improvement of calculation times (thanks to the Onelevel for example), better vision of the IOC problem, better understanding of the reliability of the results obtained and generalization of the application to human movement.

V. CONCLUSION This paper presents our holistic approach called Projected Inverse Optimal Control. This approach is complete and transversal in the framework of uncertainty related to the observation of human movement. It is to be understood that it is not a simple method of resolution. It proposes different elements, each independent of the others, to improve expression and better understand the resolution of the IOC's problem. We have shown in this paper that our approach overcomes some difficulties. The parameterization of the trajectories allows to reduce the dimensions of the problem. In the polynomial framework, we even show that there is a minimum degree depending on the state, the constraints and the number of cost functions. We also propose to decouple the conditions of the problem from its solution. The identifiability condition verifies that the problem is well posed, the singularity condition verifies an optimality condition necessary to solve the problem and the feasibility condition verifies the reliability of the results. Finally, we show that the solution of the IOC problem could be illustrated simply by projections in the space of trajectories. This vision of the resolution allowed us to explain the differences between the classical resolution methods of the literature as well as to discuss an algorithm allowing to choose a base of cost functions.

Fig. 4 :

 4 Fig. 4: Illustration of 3 singular curves: C = [τ 1 , τ 2 ] (yellow), C ′ = [ θ1 τ 1 , θ2 τ 2 ] (blue) and C ′′ = [ θ1 τ 1 , θ2 τ 2 ] (purple). The green point correspond to a uncertain trajectory. It is not obvious which cost function basis to choose in this configuration.


Joseph Laureanti

Juan Brandi

Elvis Offor

David Engel

Robert Rallo

Bojana Ginovska

Xavier Martinez

Marc Baaden

Nathan A Baker
email: nathan.baker@pnnl.gov

Visualizing biomolecular electrostatics in virtual reality with UnityMol-APBS

Keywords: electrostatics, solvation, molecular visualization, virtual reality

Virtual reality is a powerful tool with the ability to immerse a user within a completely external environment. This immersion is particularly useful when visualizing and analyzing interactions between small organic molecules, molecular inorganic complexes, and biomolecular systems such as redox proteins and enzymes. A common tool used in the biomedical community to analyze such interactions is the Adaptive Poisson-Boltzmann Solver (APBS) software, which was developed to solve the equations of continuum electrostatics for large biomolecular assemblages. Numerous applications exist for using APBS in the biomedical community including analysis of protein ligand interactions and APBS has enjoyed widespread adoption throughout the biomedical community. Currently, typical use of the full APBS toolset is completed via the command line followed by visualization using a variety of two-dimensional external molecular visualization software. This process has inherent limitations: visualization of three-dimensional objects using a two-dimensional interface masks important information within the depth component. Herein, we have developed a single application, UnityMol-APBS, that provides a dual experience where users can utilize the full range of the APBS toolset, without the use of a command line interface, by use of a simple graphical user interface (GUI) for either a standard desktop or immersive virtual reality experience.

Introduction

Understanding and predicting biomolecular processes requires a basic knowledge of the electrostatic interactions possible within a molecular system [START_REF] Fried | Electric fields and enzyme catalysis[END_REF][START_REF] Jurrus | Improvements to the apbs biomolecular solvation software suite[END_REF][START_REF] Dolinsky | PDB2PQR: expanding and upgrading automated preparation of biomolecular structures for molecular simulations[END_REF]. Many programs1 -such as PyMOL [START_REF] Schrödinger | The PyMOL molecular graphics system, version 1.8[END_REF], VMD [START_REF] Humphrey | VMD: visual molecular dynamics[END_REF], Chimera [START_REF] Pettersen | UCSF Chimera-a visualization system for exploratory research and analysis[END_REF], and UnityMol [START_REF] Doutreligne | Unitymol: Interactive scientific visualization for integrative biology[END_REF]-provide platforms for the visualization and analysis of protein systems. However, visualization of molecular systems is currently largely completed using a variety of tools that all contain a common, and inherent limitation: three-dimensional objects are observed using a two-dimensional interface. Using a two-dimensional interface to manipulate and understand three-dimensional objects removes an important depth component that is crucial to understanding spatial resolution within a three-dimensional object. Furthermore, manipulation of multiple objects in a three-dimensional space is currently performed using a cumbersome combination of mouse or track-pad and keyboard interactions. Interacting with a three-dimensional space via a two-dimensional interface is a limitation when attempting to simultaneously investigate multiple objects such as protein-protein, protein-cofactor, or enzyme-substrate interactions. In this paper, we describe a new UnityMol-APBS virtual reality (VR) platform for visualizing these threedimensional systems. We will show brief examples of how a student or researcher would use the APBS and PDB2PQR tools to generate and visualize electrostatic information and describe how this use is facilitated by the VR interface.

Molecular visualization with UnityMol-APBS

Many popular platforms exist to visualize molecular systems via two-dimensional environments; e.g., PyMOL, VMD, Chimera, and UnityMol. Several new visualization methodologies have been implemented to provide an augmented reality (AR) or VR interface [START_REF] Sommer | From virtual reality to immersive analytics in bioinformatics[END_REF][9][START_REF] O'connor | Sampling molecular conformations and dynamics in a multiuser virtual reality framework[END_REF][11][12][13]. VR interfaces-as provided by UnityMol, Chimera (see ChimeraX) [START_REF] Goddard | Ucsf chimerax: Meeting modern challenges in visualization and analysis[END_REF], Nanome [9], and Nano Simbox [START_REF] O'connor | Interactive molecular dynamics in virtual reality from quantum chemistry to drug binding: An open-source multiperson framework[END_REF]-allow students or researchers the ability to be completely immersed within a given molecular system, thereby facilitating insight into interactions between proteins and enzymes, small organic molecules, and synthetic organometallic complexes. Single-user VR software varies from viewing-only applications such as Educhem-VR [12] to tools that provide various degrees of user interaction both in a standard desktop environment and using a VR interface; e.g., UnityMol-APBS and ChimeraX [START_REF] Goddard | Ucsf chimerax: Meeting modern challenges in visualization and analysis[END_REF]. Collaborative VR software includes Nanome [9], a collaborative VR environment for drug discovery, and Nano Simbox [START_REF] O'connor | Interactive molecular dynamics in virtual reality from quantum chemistry to drug binding: An open-source multiperson framework[END_REF], a collaborative VR environment for interacting with molecular dynamics (MD) simulations. As noted by the team behind Nano Simbox, VR environments have been shown to increase the rate of scientific discovery. For example, the developers of Nano Simbox attempted to quantify the positive impact of using a virtual environment for general molecular visualizations and manipulations and found a 10-fold increase in productivity [START_REF] O'connor | Interactive molecular dynamics in virtual reality from quantum chemistry to drug binding: An open-source multiperson framework[END_REF].

UnityMol-APBS has been developed using the Unity game engine for both standard desktop (Windows, Mac OS X, and Linux) and VR (Windows-only) interfaces. A Python terminal, implemented using IronPython, has been incorporated to provide an improved user experience through additional command-line functionality. The original UnityMol software was introduced in 2013 and was developed by Marc Baaden and his research group at CNRS in the Institute of Physico-Chemical Biology. UnityMol is available free of charge and can be downloaded from a repository at SourceForge2 . The UnityMol-APBS software described in this paper is derived from UnityMol.

Molecular electrostatics with APBS

As discussed in our article from the previous Tools for Protein Science special issue [START_REF] Jurrus | Improvements to the apbs biomolecular solvation software suite[END_REF], APBS is one of several software packages for calculating electrostatic properties of biomolecular systems3 . The APBS software was developed to solve the Poisson-Boltzmann (PB) equation [START_REF] Fixman | The Poisson-Boltzmann equation and its application to polyelectrolytes[END_REF][START_REF] Grochowski | Continuum molecular electrostatics, salt effects, and counterion binding-a review of the Poisson-Boltzmann theory and its modifications[END_REF][START_REF] Lamm | The Poisson-Boltzmann equation[END_REF] for the calculation of biomolecular solvation properties and electrostatic interactions:

-∇ ⋅ ∇φ - M i c i q i e -β(q i φ+V i ) = ρ (1)
which is solved on a domain Ω ⊆ R 3 with φ specified on the boundary of that domain. In Eq. 1,

φ ∶ Ω ↦ R is the electrostatic potential, ∶ Ω ↦ [ε u , ε v ] is a dielectric coefficient function that ranges between solute ε u > 0 and solvent ε v > 0 dielectric values, and ρ ∶ Ω ↦ R is a charge distribution function.
For each mobile ion species, i = 1, . . . , M , q i ∈ R is the charge, c i > 0 is the concentration, and V i ∶ Ω ↦ R is the steric ion-solute interaction potential. Finally, β = (kT ) -1 > 0 is the inverse thermal energy where k is the Boltzmann constant and T is the temperature. The main input to APBS is a PQR file that specifies the position in Ω, charge values for ρ, and radii (used to construct and V i ) for each atom. The PDB2PQR software was developed to expedite the preparation of input files for analysis by APBS [START_REF] Dolinsky | PDB2PQR: expanding and upgrading automated preparation of biomolecular structures for molecular simulations[END_REF][START_REF] Dolinsky | PDB2PQR: an automated pipeline for the setup of Poisson-Boltzmann electrostatics calculations[END_REF]. Functions of the PDB2PQR software includes reading a PDB file, converting to the PQR format, repairing missing heavy atoms, optimizing titration states, adding missing hydrogen atoms, assigning radius and charge parameters, and automatically preparing subsequent APBS electrostatic calculations.

APBS and PDB2PQR can be used via command-line operations, through popular molecular visualization software (Chimera, PyMOL, VMD), or via a web-based interface. As described in this article, APBS and PDB2PQR tools have now also been incorporated into UnityMol-APBS with three primary goals:

• offer an intuitive visual interface to APBS functionality without the requirement of compiling,

• enable users to employ one application to calculate and immediately visualize and compare multiple results in a completely immersive virtual reality experience, and

• provide the terminal commands, that were generated in UnityMol-APBS, to PDB2PQR and APBS as input, in a text file for the purposes of reproducibility of results.

Installation of the integrated UnityMol-APBS and APBS environment is straightforward: both executables should be downloaded from https://github.com/Electrostatics/VR and installed on the user system using default values4 .

Methods and results

The full APBS toolset, including PDB2PQR, was integrated with UnityMol-APBS to allow users access the functionality of this software without the requirement of manually compiling the tools The main menu for interactions with APBS and PDB2PQR. Blue boxes show the buttons to launch PDB2PQR and APBS executables, green boxes show the location of the options used for producing the images in Figures 2 and3, and the purple boxes highlight the two input fields required to use custom force fields and custom residue names.

and using terminal commands. UnityMol was specifically chosen for integration with APBS due to its high-quality visual representations, ease of interactions with multiple molecular objects, powerful tools for visual manipulations, and the ability to switch between a desktop and a virtual environment. All APBS tools are directly accessed from within UnityMol-APBS (see Figure 1) and can be used from within the desktop application or a virtual environment to allow immediate comparison of many structural properties including titration states, electrostatic surface potentials, electrostatic field lines, hydrogen bonding interactions, and salt bridge formation.

The examples below illustrate how to visualize the properties of proteins with natural and unnatural residues. Values and methodologies used in the examples are intended as a demonstration and are not necessarily the recommended values for all applications. Please refer to the PDB2PQR and APBS documentation when using APBS for your own applications. The HTC Vive VR system was used for the examples shown in this paper; however, several VR systems are supported including: Oculus CV1, DK2, and Rift S. Both the desktop and VR environments of UnityMol-APBS save all user commands as a Python script that can be loaded later to reproduce or reuse the visualization workflow. This greatly reduces the time required to begin working from where a researcher left off, or when a student is loading a scene that an instructor has previously prepared.

Analyzing acetylcholinesterase (AChE) electrostatic surface potentials and field lines

Torpedo californica AChE is an enzyme that operates incredibly fast with high selectivity, in part, due to its electrostatic properties [START_REF] Quinn | Acetylcholinesterase: Enzyme structure, reaction dynamics, and virtual transition states[END_REF]. Catalysis occurs at an active site buried within the enzyme. Electric fields help guide positively charged substrates from the enzyme exterior through a narrow channel to the buried active site as shown in the UnityMol-APBS visualization in Figure 2. As such, this enzyme is a good example for exploring electrostatic surface potentials and electrostatic field lines in an immersive and virtual experience. The instructions provided here for the AChE analysis can be used as a template for other uses of UnityMol-APBS. The Appendix includes step-by-step instructions for this example (Example A.1) as well as the output Python code (Example A.2), which can be used by a user to instantly load our results. Structural information for Torpedo californica AChE was obtained from PDB entry 5EI5, via the fetch command in UnityMol-APBS. Once loaded, the default force field parameters for AM-BER [START_REF] Case | The AMBER biomolecular simulation programs[END_REF] were chosen and amino acid protonation states were estimated using PROPKA at pH 7.0. PROPKA is a heuristic method for computing pK a values by incorporating effects due to desolvation, hydrogen bonding, and charge-charge interactions [START_REF] Sondergaard | Improved treatment of ligands and coupling effects in empirical calculation and rationalization of pK a values[END_REF]. Figure 1 highlights the options used with PDB2PQR: -chain, -apbs-input, -summary, -drop-water, -salt, and -hbond, to retain the chain information, prepare an input file for APBS, print a summary, remove water molecules, write salt bridge interactions, and write hydrogen bonding interactions, respectively. UnityMol-APBS creates and saves all files in the destination directory C:\APBS PDB2PQR\OutputFiles. Upon ) with bound alkylene-linked bis-tacrine. (A) Electrostatic field lines and protein secondary structure shown with alpha helices (yellow), beta sheets (blue), and random coils (white). Residues Tyr70, Trp84, Trp279, and Phe330 are shown interacting with alkylenelinked bis-tacrine via hydrogen bonding and π-π stacking interactions. The red oval highlights the potential gradient. (B) AChE surface model with field lines and mapped electrostatic surface potentials shown with neutral, negative, and positive charges in white, red, and blue, respectively. Field lines are calculated from a gradient (value of 0.2) and depicted with the starting points in red and the ending points in blue. The orientation is the same in Figures A andB, where the alkylene-linked bis-tacrine can be seen occupying the catalytic gorge. The white circle highlights the potential gradient exiting the catalytic gorge.

running the PDB2PQR tools, the newly created .PQR file is immediately loaded for the user. Hydrogen bonding and salt bridge information are written as .hbonds and .salt formatted files; this information can be read and visualized directly from within UnityMol-APBS as independent selections. This allows users to color the interactions as desired.

Default parameters for APBS Poisson-Boltzmann calculations were selected through the Unity-Mol-APBS interface. APBS writes electrostatic potential output in OpenDX format that is directly loaded by UnityMol-APBS after completing the calculation. The potential gradient within the OpenDX file allows UnityMol-APBS to color the molecular surface based on electrostatic potential and show the associated electrostatic field lines, Figure 2. Highlighted in Figure 2 is the strong electric field originating in the AChE active site channel. Due to the confinement of the catalytic gorge, this channel is extremely challenging to explore with traditional molecular visualization methods. Additionally, manually moving a ligand through the catalytic gorge to analyze the residues that could interact is not a task that can be completed easily in 2D visualization environments. Completing this task in VR is much easier as the depth component is not obscured due to the inherent limitations associated with viewing a three-dimensional space on a two-dimensional screen. Furthermore, the virtual environment more easily allows the user to simultaneously manipulate both the enzyme and the substrate to provide optimized viewing angles. Animated electrostatic field lines in UnityMol-APBS represent the electrostatic potential gradient and the following properties can be manipulated: movement speed of the lines, width of lines, length of the line, and the gradient used to calculate the lines. An example of the electrostatic surface potential and electrostatic field lines of Torpedo californica AChE is shown in the supplementary information video, which can be accessed at https://github.com/Electrostatics/VR. Also, included in the video is an example of a user manipulating two objects (substrate and AChE) simultaneously in a virtual space manually traversing the substrate through the catalytic gorge of AChE.

Using custom force fields to analyze an artificial enzyme: Rhodium-containing Lactococcal Multidrug Resistance Regulator (Rh-LmrR*)

The electrostatic field lines an electrostatic surface potentials of an artificial metalloenzyme capable of hydrogenating CO 2 to a liquid fuel, formate, is shown in Figure 3 [START_REF] Laureanti | Protein scaffold activates catalytic co2 hydrogenation by a rhodium bis(diphosphine) complex[END_REF]. This system, Rh-LmrR (where LmrR is the Lactococcal multidrug resistance regulator and Rh refers to an organometallic rhodium-(bis)diphosphine complex), was chosen as the second example to illustrate the use of custom force fields with UnityMol-APBS. The artificial metalloenzyme contains a covalently anchored Rh(bis)diphosphine cofactor that cannot be analyzed using standard force field parameters due to the non-standard rhodium and phosphorus atoms, as well as the associated ligand framework.

In this example, we show how to input externally generated force field data, analyze hydrogen bonding interactions, identify salt bridge locations, calculate the electrostatic potential at multiple atoms using the APBS multivalue tool, as well as visualize the electrostatic surface potential and electrostatic field lines. Our methods to prepare the custom force fields, using tools outside of UnityMol-APBS, are briefly described here. Structural data for the artificial metalloenzyme was first obtained from crystallographic data (PDB ID 6DO0) with further minimization via density functional theory (DFT) calculations and conformational sampling using MD simulations in GROMACS [24]. The resulting PDB file for Rh-LmrR was loaded into UnityMol-APBS and PDB2PQR was used to prepare the structure for APBS calculations as described in the AChE example above. The parameters for the Rh-(bis)diphosphine complex do not exist in any of the included standard libraries, therefore custom force field parameters for Rh-LmrR were explicitly input to UnityMol-APBS. Calculations for DFT were performed in NWChem 6.5 [START_REF] Valiev | NWChem: A comprehensive and scalable open-source solution for large scale molecular simulations[END_REF] using the B3LYP functional [START_REF] Becke | Densityfunctional thermochemistry. III. The role of exact exchange[END_REF][START_REF] Stephens | Ab initio calculation of vibrational absorption and circular dichroism spectra using density functional force fields[END_REF] and 6-31G* basis set [START_REF] Petersson | A complete basis set model chemistry. I. The total energies of closedshell atoms and hydrides of the firstrow elements[END_REF][START_REF] Petersson | A complete basis set model chemistry. II. Openshell systems and the total energies of the firstrow atoms[END_REF] for all atoms except for Rh, which used a large-core correlation-consistent pseudopotential basis set [START_REF] Martin | Correlation consistent valence basis sets for use the Stuttgart-Dresden-Bonn relativistic effective core potentials: The atoms Ga-Kr and In-Xe[END_REF]. The geometries were optimized in the gas phase and verified as minima using Hessian calculation in a rigid-rotor, harmonic approximation. For the charge calculation, the complex was capped with methyl groups, constrained to an overall neutral charge, and all equivalent atoms were constrained to equal values. Charges for atoms other than Rh and P were calculated using the RESP procedure [START_REF] Cornell | Application of resp charges to calculate conformational energies, hydrogen bond energies, and free energies of solvation[END_REF]. The atomic point charges of the Rh center and phosphorus atoms were obtained using the Blöchl scheme [START_REF] Blöchl | Electrostatic decoupling of periodic images of plane-wave-expanded densities and derived atomic point charges[END_REF], since the RESP scheme resulted in unphysical values. Lennard-Jones parameters for each atom were taken from the GAFF force field [START_REF] Wang | Development and testing of a general amber force field[END_REF], except for Rh which used UFF [START_REF] Rappe | UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations[END_REF].

UnityMol-APBS allows custom force fields and residue names to be used by duplicating the base force field files and amending new custom force field and residue name data as described in the APBS-PDB2PQR documentation. Once these files are prepared, the data can be loaded via the PDB2PQR menu in UnityMol-APBS. The location of the custom force field is provided by use of the "Find FF" button, and the custom residue names are found by use of the "Find FF names" button, shown highlighted in purple in Figure 1. When the PDB2PQR executable is launched via UnityMol-APBS (highlighted blue in Figure 1), the paths to the custom files will be automatically loaded, overriding the force field information chosen at the top of the GUI. After the PDB2PQR process finishes, the newly generated .pqr file is automatically loaded into UnityMol-APBS as a new selection. The APBS executable (highlighted in blue, Figure 1) can then be launched to solve the Poisson-Boltzmann equations for continuum electrostatics. The newly generated OpenDX file will be automatically loaded into UnityMol-APBS, thus facilitating visualization of the electrostatic surface potential and electrostatic field lines.

Please see attached image. Using a virtual environment to complete this task is simpler than a two-dimensional interface. Dihedral angles can also be shown directly using the UnityMol-APBS interface. Distance measurement line and angle measurement line are shown as blue lines. Grey spheres are used to highlight the selected atoms.

As an example of how to evaluate specific atomic potential values from within UnityMol-APBS, we used the immersive virtual reality interface to compare the electrostatic potential of the four phosphorus atoms and the Rh metal center of the organometallic complex alone in solution to the Rh complex immobilized within the protein scaffold. In the original publication [START_REF] Laureanti | Protein scaffold activates catalytic co2 hydrogenation by a rhodium bis(diphosphine) complex[END_REF], it was found that the Rh-complex alone in solution was not a competent catalyst for CO 2 hydrogenation. However, catalytic activity was observed once the Rh complex was immobilized in the protein scaffold. The protein scaffold imposes conformational restrictions and as this example implies, also perturbs the electrostatic potential around the Rh-P environment. This shift in electrostatic potential was observed using the APBS multivalue tool via UnityMol-APBS. To complete this task without imposing influences from the grid, we used the following workflow to accurately calculate potentials while removing "self-energy" artifacts from the grid. First, two initial .dx files containing the calculated electrostatic potential map were created from the respective (APBS) input files for the Rh-complex and the artificial enzyme using a protein dielectric (pdie) of 2 and a solvent dielectric (sdie) of 78.4. Second, an additional .dx file was created for each system with a homogeneous dielectric of 2 (pdie=sdie=2). Third, the data from each .dx file was then used to generate potential values from the multivalue tool by selecting the desired atoms (one Rh atom and four P atoms) of the active site. Fourth, the reaction field potential values were calculated by subtracting the homogeneous-dielectric values (second step) from the inhomogeneous-dielectric values (first step) to remove self-energy grid artifacts in the calculation. Fifth, the APBS coulomb tool was used to analytically calculate the electrostatic potential at each atom, assuming a homogeneous environment of dielectric 2. The analytically calculated Coulomb potential results were added to the reaction field potential results to give the final potential value at each atom site. From these results, we observed that all four phosphorus atoms were perturbed upon immobilization with an average shift of 165 ± 15 mV relative to the Rh-complex alone in solution. The Rh center shifted -170 mV upon immobilization within the protein environment.

Using VR to explore the active site of an artificial enzyme in an immersive environment is extremely helpful when trying to understand the electrostatic environment, spatial constraints, local flexibility, and potentially beneficial or inhibiting interactions of neighboring residues. A full immersion experience allows the user to observe the protein surface from within the protein scaffold, a feature that is drastically hindered when using a standard two-dimensional interface. As illustrated in Figure 4, the VR environment is particularly well-suited for measuring distances, angles, and electrostatics potentials at specific sites in the protein (e.g., the Rh metal center). The APBS multivalue tool can be used to identify the electrostatic potential at specific atoms selected in UnityMol-APBS and is described in further detail in the Appendix (Example A.3).

Conclusions

In this newly coupled version of UnityMol-APBS, we have provided a seamless desktop/virtual interface to provide a robust method of interaction with molecular systems. The virtual world allows a user to easily manipulate and compare multiple three-dimensional objects in real time, a task that is very much hindered using only a standard desktop interaction. Additionally, the virtual interface allows users to quickly prepare input files for the PDB2PQR and APBS toolset by simply reaching out and grabbing desired atoms, which are then automatically prepared by the UnityMol-APBS software and delivered to the APBS tools via command line interactions. This dramatically reduces the time required, since the user no longer has to sift through multiple text files to prepare custom input files for the PDB2PQR and APBS toolsets. Furthermore, the inherent need to interact with the PDB2PQR and APBS tools through direct command line interactions has been removed. All command line inputs are setup using a helpful GUI and delivered directly from UnityMol-APBS. This greatly increases the availability of the PDB2PQR and APBS toolsets to undergraduate students as well as anyone entering the field of molecular visualization, as users are not required to compile the software before use, nor is there an inherent need to understand command line interactions.

APBS has been in public use for nearly 20 years, thanks to ongoing support by the National Institutes of Health to maintain and update the codes, work with the user base to ensure access and execution of the codes, and improve performance and incorporate new features based on new algorithms and user feedback. Any software package still in use over that span of time requires continuous updating and refinement to leverage emerging software technologies and improvements in measurement technology that affect the input data. Therefore, APBS has been redesigned, as illustrated in Figure 5, to ensure the continued availability of this free and scalable software package for biomolecular electrostatics, solvation, and structure assessment/preparation. The goal of the ongoing software re-engineering effort is to unify the software development into a single code-base that can be easily deployed on desktops, servers and cloud infrastructures. To achieve this goal, the software has been refactored as a microservices architecture to improve its extensibility and scalability. The redeveloped software is currently under testing. Access to this new version of APBS will be available through http://www.poissonboltzmann.org; the underlying developmental code is available at https://github.com/Electrostatics/apbs-rest.

Example A.2 (Output Python code from UnityMol-APBS for AChE VR visualization in Section 2.1). load("D:/5EI5_May/5ei5.pdb", True, False) load("D:/5EI5_May/5ei5X.pqr", True, False) loadDXmap("5ei5X", "D:/5EI5_May/5ei5X.dx") showSelection("all(5ei5X)", "fl") showSelection("all(5ei5X)", "s") colorByCharge("all(5ei5X)", "s") #UnityMol Specific functions - ---------- 9. Select the "+" button on the "all(RhLmrRX)" selection tab and select "surface" 10. Select "color by charge" under the surface selection tab 11. Select the "+" button on the "all(RhLmrRX)" selection tab and select "field lines" 12. Select the "New Selection" button under the "Utils" tab 13. Click "Selection mode" twice to change to "Atom" 14. Select the four phosphorus atoms and the Rh atom, the new selection should contain five atoms 15. With the new selection highlighted, click the "multivalue" tool button, under the APBS tools tab, to evaluate the electrostatic potential for these five atoms.

Two important text files are created: 1) containing the original data for these atoms from the .pdb file and 2) the output from multivalue. The order of the lines in file 1 are equivalent to that of file 2. This is helpful when keeping track of multiple atoms evaluated with multivalue.

Figure 1 :

 1 Figure 1: UnityMol-APBS user interface for PDB2PQR and APBS. (A) The main UnityMol-APBS menu; orange box highlights the two buttons used to open the APBS and PDB2PQR tools. (B)The main menu for interactions with APBS and PDB2PQR. Blue boxes show the buttons to launch PDB2PQR and APBS executables, green boxes show the location of the options used for producing the images in Figures2 and 3, and the purple boxes highlight the two input fields required to use custom force fields and custom residue names.

Figure 2 :

 2 Figure2: Electrostatic surface potential and field lines of Torpedo californica AChE (PDB ID 5EI5) with bound alkylene-linked bis-tacrine. (A) Electrostatic field lines and protein secondary structure shown with alpha helices (yellow), beta sheets (blue), and random coils (white). Residues Tyr70, Trp84, Trp279, and Phe330 are shown interacting with alkylenelinked bis-tacrine via hydrogen bonding and π-π stacking interactions. The red oval highlights the potential gradient. (B) AChE surface model with field lines and mapped electrostatic surface potentials shown with neutral, negative, and positive charges in white, red, and blue, respectively. Field lines are calculated from a gradient (value of 0.2) and depicted with the starting points in red and the ending points in blue. The orientation is the same in FiguresA and B, where the alkylene-linked bis-tacrine can be seen occupying the catalytic gorge. The white circle highlights the potential gradient exiting the catalytic gorge.

Figure 3 :

 3 Figure 3: Single conformation from a MD simulation of Rh-LmrR* showing (A) the protein secondary structure and the immobilized Rh-complex. (B) electrostatic surface potential and field lines. Color scheme is identical to Figure 2.

Figure 4 :

 4 Figure 4: Example of using UnityMol-APBS to explore atomic distances and bond angles.Using a virtual environment to complete this task is simpler than a two-dimensional interface. Dihedral angles can also be shown directly using the UnityMol-APBS interface. Distance measurement line and angle measurement line are shown as blue lines. Grey spheres are used to highlight the selected atoms.

Figure 5 :

 5 Figure 5: Diagram for redesigned and refactored APBS software to enable modular growth.

Figure 1 .Figure 2 .Figure 3 .Figure 4 .Figure 5 .

 12345 Figure 1. UnityMol-APBS user interface for PDB2PQR and APBS. (A) The main UnityMol-APBS menu; orange box highlights the two buttons used to open the APBS and PDB2PQR tools. (B)The main menu for interactions with APBS and PDB2PQR. Blue boxes show the buttons to launch PDB2PQR and APBS executables, green boxes show the location of the options used for producing the images in Figures2 and 3, and the purple boxes highlight the two input fields required to use custom force fields and custom residue names.179x161mm (150 x 150 DPI)

See https://en.wikipedia.org/wiki/List_of_molecular_graphics_systems for a partial list of molecular visualization programs.

UnityMol can be downloaded at https://sourceforge.net/projects/unitymol.

The APBS software can be downloaded or used via the web at http://www.poissonboltzmann.org.

UnityMol-APBS assumes a default installation of APBS at C: APBS PDB2PQR on Windows systems.

Acknowledgments

The authors gratefully acknowledge NIH grant GM069702 for support of APBS and PDB2PQR. The development of the VR methods was supported by the US Department of Energy (DOE), Office of Science, Office of Basic Energy Sciences (BES), Division of Chemical Sciences, Geosciences, & Biosciences. MB acknowledges support by the "Initiative d'Excellence" program from the French State (Grant "DYNAMO", ANR-11-LABX-0011, and "CACSICE", ANR-11-EQPX-0008). PNNL is operated by Battelle for the U.S. DOE under contract DE-AC05-76RL01830.

A Appendix

The following examples provide step-by-step directions for the systems described in the main text. The source files are also available for download as Supporting Materials for this article at https: //github.com/Electrostatics/VR. Select "all(5EI5)" and run PDB2PQR 6. 5ei5X.pqr is written to a file and is immediately loaded for the user.

7. Select "all(5EI5)" and run APBS 8. 5ei5X.dx is created and loaded into the selection "all(5EI5X)" automatically 9. Select the "+" button on the "all(5EI5X)" selection tab, then select "surface" 10. Select "color by charge" 11. Select the "+" button on the "all(5EI5X)" selection tab, then select "field lines"


30 > BIOFUTUR 386 • MAI 2017 de la cryoablation de tumeur rénale ou l'embolisation d'anévrisme, constituent, en effet, une alternative à la chirurgie. D'autres sont un complément thérapeutique d'un traitement médical parfois inefficace, la consolidation d'une vertèbre tassée hyperalgique par cimentoplastie par exemple. Enfin, certaines se présentent comme la seule thérapie envisageable, comme la thrombectomie mécanique dans le traitement de l'accident vasculaire cérébral (AVC) ou encore l'embolisation d'hémostase. * 1 Spires métalliques utilisées comme matériel d'embolisation (occlusion vasculaire) * 2 Exposition à une température donnée pendant un temps donné

RADIOLOGIE INTERVEN-TIONNELLE VASCULAIRE : L'EXEMPLE DE L'AVC

La radiologie interventionnelle vasculaire utilise les vaisseaux sanguins pour faire cheminer du matériel d'un point de ponction cutané -le plus souvent fémoralvers une cible accessible par lesdits vaisseaux. Ce chemin naturel permet de déposer du matériel dans un vaisseau endommagé (prothèse dans une sténose, coils* 1 dans un anévrisme), de réaliser une chimiothérapie ciblée (chimio-embolisation hépatique), mais aussi de désobstruer une artère occluse.

La thrombectomie mécanique a récemment été validée dans la prise L a radiologie interventionnelle a pour finalité le traitement de diverses pathologies (tumeurs, lésions vasculaires) par des techniques dites « minimalement invasives » guidées par l'imagerie, en y accédant par les vaisseaux ou par la peau. Ses débuts sont attribués à Charles Dotter qui, à la tête du département de radiologie vasculaire et interventionnelle de l'Université d'Oregon, eut l'idée de franchir le cap de la fonction purement exploratrice de l'imagerie invasive. Dès 1963, en effet, il proposa d'emprunter la voie vasculaire, déjà utilisée pour l'imagerie diagnostique invasive par insertion de cathéters, pour traiter des lésions vasculaires.

Le développement de l'imagerie médicale (imagerie par rayons X, échographie et IRM) ouvrit, dès les années 1980, la porte à de nombreux champs d'applications thérapeutiques. Les procédures minimalement invasives connaissent un succès grandissant depuis la fin des années 1990 puisqu'elles s'accompagnent d'une réduction de la durée d'hospitalisation tout en permettant de traiter des patients pour lesquels la chirurgie est contreindiquée, avec, en corolaire, une diminution de la morbidité. Certaines de ces procédures, à l'image en charge de l'AVC ischémique (1,2). La privation d'oxygène des neurones consécutive à l'occlusion brutale d'une artère peut engendrer la nécrose d'une partie du cerveau en quelques minutes. Par opposition à cette nécrose irréversible, les neurones peuvent rester dans un état de sidération réversible grâce à la restauration du flux. C'est là qu'intervient le radiologue, qui va retirer le thrombus dans l'artère le plus tôt possible, idéalement dans les 4 à 5 heures suivant l'AVC. La technique consiste à déployer un stent dans le thrombus qui, en s'agrippant à ses mailles, va pouvoir être retiré dans le cathéter, en même temps que le stent. Plusieurs études randomisées ont montré que la restauration du flux s'accompagne d'une meilleure récupération neurologique comparativement à l'injection seule, en intraveineux, d'un agent thrombolytique (3). Plusieurs aiguilles peuvent être utilisées simultanément pour couvrir des volumes de plusieurs dizaines de cm 3 , avec l'avantage de bien maîtriser visuellement les limites de la congélation au scanner ou à l'IRM. Les indications sont principalement les tumeurs du rein, du poumon et certaines lésions rares des tissus mous (tumeurs desmoïdes). Quelle que soit la technique utilisée, l'imagerie médicale est essentielle pour s'assurer du succès de la thérapie. Les procédures d'ablation thermique guidées par imagerie ont été standardisées et sont découpées en plusieurs étapes clés. En premier lieu, la planification, dont l'objectif est de repérer le nombre et la localisation des cibles, d'anticiper le nombre d'applicateurs nécessaires et de déterminer la meilleure balistique. Ensuite, le ciblage, consistant à amener chaque applicateur vers la zone désirée en utilisant l'imagerie en temps réel. Enfin, le monitoring en temps réel pendant l'ablation, dont l'objectif est de s'assurer qu'aucun dommage imprévu ne vient toucher les tissus environnants, tout en contrôlant qu'une dose thermique suffisante et nécessaire soit appliquée à la lésion.

RADIOLOGIE INTERVEN-TIONNELLE PERCUTANÉE : L'EXEMPLE DE L'ABLATION THERMIQUE

VERS UNE THÉRAPIE NON INVASIVE ?

Une nouvelle technique d'ablation, qui diffère des autres dans la mesure où elle est absolument non invasive, illustre une piste de développement de la radiologie interventionnelle dans le domaine de l'oncologie. La thérapie par ultrasons focalisés de haute intensité (HIFU : High Intensity Focused Ultrasound) repose sur l'émission de faisceaux ultrasonores focalisés sur la zone à traiter et sur les effets mécaniques et/ou thermiques des ultrasons sur les tissus. Cela permet, par exemple, d'élever la température localement autour d'un foyer afin d'y induire une lésion thermique. Cette thérapie est déjà utilisée depuis plusieurs années pour le traitement des fibromes utérins ou du cancer de la prostate (9-13). La communauté de la radiologie interventionnelle est en passe d'adopter cette technique prometteuse, présentée par certains comme de la « chirurgie non invasive ». L'analogie se fait rapidement avec la radiothérapie stéréotaxique, version également non invasive, plus ancienne et ionisante de destruction tumorale, qui partage la phase de planification avec l'imagerie interventionnelle.

Les défis de la radiologie interventionnelle sont donc à la fois technologiques, économiques et scientifiques. La radiologie interventionnelle est, en effet, intime-

Cryoablation réalisée sous IRM

Procédure réalisée dans le service de radiologie interventionnelle des Hôpitaux universitaires de Strasbourg. Ces procédures nécessitent une expertise pointue au sein même de la communauté de la radiologie interventionnelle, un matériel adapté compatible avec l'IRM, et des séquences IRM dites « interventionnelles », particulièrement rapides et interactives.

Thérapie par HIFU

La thérapie par HIFU permet de cibler une zone directement en profondeur, de manière non invasive. Le contrôle par IRM permet notamment d'évaluer la thérapie via une mesure de température en temps réel.

  Les interventions percutanées concernent à la fois les gestes diagnostiques -prélèvement d'une lésion pour en établir les causes -et thérapeutiques, avec des champs d'application très variés : infiltrations de corticoïdes, sclérothérapie d'une malformation vasculaire, destructions thermiques d'une lésion… Les ablations thermiques visent à détruire un tissu cible -une tumeur, par exemple -par hyper ou hypothermie, en induisant une nécrose et/ou une apoptose cellulaire. La grandeur « dose thermique »* 2 a été proposée dans les années 1980 en tant que métrique permettant de quantifier cette mort cellulaire. Par exemple, sur la majorité des tissus, une exposition de plusieurs dizaines de minutes à 43 °C divise par 10 le nombre de cellules vivantes. Cette relation est fortement non linéaire puisqu'il ne suffit que de quelques minutes à 60 °C pour obtenir le même effet (4). En radiologie interventionnelle, l'ablation est réalisée à l'aide d'un applicateur, une aiguille intégrant le dispositif permettant d'élever ou d'abaisser la température du tissu. Cela permet donc de réaliser des ablations ciblées, potentiellement sur toute partie du corps, de manière minimalement invasive. Parmi les techniques d'ablation thermique les plus courantes, citons : • l'ablation par radiofréquences : l'une des plus courantes et des plus anciennes, elle repose sur l'échauffement du tissu résultant du passage d'un courant généré par des électrodes (5). La taille des lésions thermiques est typiquement de 2 à 4 cm 3 . Les lésions tumorales focalisées du foie ou du poumon sont de bonnes indications ; • l'ablation par micro-ondes : plus récente, cette technique repose sur l'absorption par les tissus de l'énergie électromagnétique générée par l'applicateur émettant des micro-ondes entre 900 mégahertz et 2,5 gigahertz (6). Une zone assez volumineuse peut être chauffée (> 5 cm 3 ). Les indications sont sensiblement identiques à celles de la radiofréquence ; • l'ablation par laser : l'applicateur se présente sous la forme d'une fibre optique insérée au centre de l'aiguille creuse permettant sa mise en place (7). La puissance utilisée est de l'ordre de quelques watts pendant 1 à 20 minutes. Une application courante est l'ablation de lésions osseuses de petit volume, comme l'ostéome ostéoïde, lésion tumorale bénigne parfois très douloureuse ; • la cryothérapie : c'est, ici, le froid, dont l'effet destructeur est connu depuis le XIX e siècle, qui est utilisé localement dans la zone à traiter. Les traitements modernes par cryothérapie reposent sur la détente d'un gaz comprimé (argon), qui, passant à pression atmosphérique, se refroidit par effet Joule-Thomson (8). Une aiguille de cryoablation est composée de plusieurs chambres permettant la détente du gaz.

  ment liée aux progrès technologiques, aussi bien ceux de l'imagerie médicale que ceux du matériel thérapeutique. Son adoption par la plupart des spécialités médicales et chirurgicales, l'intégration et la mutualisation des salles d'imagerie dans les blocs opératoires montrent qu'elle fait partie intégrante des nouvelles modalités thérapeutiques. Comme toute technologie innovante, son développement -coûteux -passe par une évaluation clinique rigoureuse, ce qui a, par exemple, été réalisé pour le traitement de l'AVC ou encore celui des tassements vertébraux. I Les actes interventionnels contrôlés par imagerie sont en pleine expansion. Cette augmentation est intimement liée aux progrès technologiques dans des champs disciplinaires allant des systèmes d'imagerie et du traitement d'images aux dispositifs d'intervention, tels que les cathéters et les robots. Smith WS et al. (2008) Stroke 39, 1205-12 (2) Jovin TGet al. (2015) N Engl J Med 372, 2296-306 (3)Saver JL et al. (2015) NEngl J Med 372, 2285-95 (4) Sapareto SA, Dewey WC (1984) Int J Radiat Oncol 10, 787-800 (5) McGahan JP et al. (1990) Invest Radiol 25, 267-70 (6) Dupuy DE (2009) Radiology 251, 617-8 (7)Nikfarjam M, Christophi C (2003) Br J Surg 90, 1033-47 (8)Buy X et al. (2013) Am J Roentgenol 201, 1353-61 (9) Goldberg SN et al. (2005) Radiology 235, 728-39 (10) Hesley GK et al. (2013) Cardiovasc Intervent Radiol 36, 5-13 (11)Souchon R et al.(2003) Ultrasound Med Biol29, 1007-15 (12) Sundaram KM et al. (2017) Semin Interv Radiol 34, 187-200 (13)Woodrum DA et al.(2016) Abdom Radiol NY 41, 877-88 Thrombectomie mécanique d'un AVC sylvien droit IRM (A. séquence diffusion ; B. séquence Time of Flight) montrant une souffrance ischémique dans le territoire sylvien profond droit (étoile) et une occlusion de la naissance de l'artère sylvienne droite. C. Artériographie cérébrale montrant l'occlusion sylvienne droite (flèche). D. Déploiement d'un stent dans la thrombose (flèche) avancé dans un microcathéter (flèche en pointillés). Le thrombus est retiré (E) et le contrôle montre une ouverture complète de l'artère sylvienne (F).

© ICUBE/HUS STRASBOURG/CH ANNECY © ICUBE/HUS STRASBOURG © ICUBE, STRASBOURG