Benjamin Michaud
email: benjamin.michaud@umontreal.ca

François Bailly
email: francois.bailly@inria.fr

Eve Charbonneau
email: eve.charbonneau.1@umontreal.ca

Amedeo Ceglia
email: amedeo.ceglia@umontreal.ca

Lea Sanchez
email: lea.sanchez@umontreal.ca

Mickael Begon
email: mickael.begon@umontreal.ca

Franc ¸ois Bailly

L Éa Sanchez

Bioptim, a Python framework for Musculoskeletal Optimal Control in Biomechanics

Keywords: Biomechanics, Musculoskeletal simulation, Optimal control, Optimization, Software

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

I. INTRODUCTION

Biomechanics researchers rely on numerical simulations of motion to gain understanding on a variety of scientific topics such as the physiological causes of movement disorders and their consequences on health [?], the estimation of nonmeasurable physiological quantities (e.g., muscle forces [?]) and the optimality of human movement [?]. The musculoskeletal models used in these simulations generally have a large number of degrees of freedom and they are governed by several ordinary differential equations (ODEs) which mainly describe multibody and muscle activation dynamics. The complexity of these systems has led scientists to formulate their simulations as optimal control problems (OCP), relying on efficient non-linear optimization software to find trajectories that fulfill a desired task while enforcing the system dynamics and minimizing a cost (e.g. motion duration, energy expenditure, matching experimental data, etc.). Up to very recently, there was no off-the-shelf software available to the community to quickly formulate and solve such musculoskeletal OCPs [?]. Consequently, researchers had to develop their own solutions, with little or no dissemination to the community, limiting synergies between researchers.

As a result, many approaches coexist to formulate and solve OCPs in the biomechanical literature. The formulation, also called transcription, consists in turning a continuous trajectory optimization problem into a generic discrete non-linear program (NLP) that is solved using a dedicated algorithm. The main family of so-called direct transcription methods comes from numerical optimal control. They consist in straightforwardly choosing the state and/or the control as optimization variables at a given number of points along the trajectory and they rely on the integration of the system dynamics between these points.

For instance, the direct collocation method has shown its efficiency in several studies investigating human motion [?],

[?], [?], [?]. It consists in approximating the integration of the system dynamics using polynomials that describe the state and control trajectories. Its main features are that it leads to very sparse NLPs, that knowledge about the state trajectory can be used in the initialization, and that it handles unstable systems well. Its major disadvantage is that adaptive integration error control implies regridding the whole problem and thus changes the NLP dimensions [?]. Direct multiple shooting is another direct method that was also applied with success in a lot of biomechanics [?], [?], [?], [?] and robotics [?], [?], [?] studies. Its features are mostly the same as for direct collocation in addition to combining absolute local error control with fixed NLP dimensions, as it relies on possibly adaptive ODE solvers to integrate the system dynamics. It leads to smaller yet less sparse NLPs than direct collocation, making the preferential use of one method over the other case-dependent. This leads to a lack of consensus in the community on which method is the most appropriate for solving biomechanical OCPs. Besides direct methods, other choices can be made, as in [?], [?], where the optimization variables are instants at which a switch in the motor strategy occurs, using polynomials function (4th, 5th order) in-between, or in [?], [?], where the optimization variables are the coefficients of fourth order polynomial approximations of the states, with linking conditions to enforce the continuity of the controls. These last approaches are less generic than the direct methods as they either require a prior knowledge about the state and control trajectories. Most of the time, when investigating complex biomechanics issues, we do not have this information.

Concerning the non-linear solver, a variety of software exist and have been used to solve musculoskeletal NLPs. They can use different heuristics: interior point methods (Ipopt, [?]) or sequential quadratic programming (SNOPT [?], Acados [?]), but they are all gradient based. Therefore, derivatives of the NLP cost function and constraints are required to perform optimization. These derivatives can be obtained by finite differences (often implemented but inaccurate thus compromising convergence) or computed exactly using automatic differentiation (requiring to write all dependencies of the software in symbolic variables), using, e.g., CasADi [?].

To promote the use of musculoskeletal optimal control in biomechanics research, we identified a strong need for a dedicated tool, as shown by the recently launched SCONE (2018) [?] and OpenSim Moco (2020) [?]. The biomechanics community being mainly composed of software users, such a tool should requests a flexible user interface written in a widely used high-level and if possible open-source language (e.g. Python) with a low-level core (e.g. C++) for efficiency. To develop such a software, four interrelated components are essential in our opinion: i) a musculoskeletal modeling software, with a visualization module (multibody kinematics and dynamics, muscle dynamics, etc.), ii) a method for automatic differentiation, iii) a discretization approach, and iv) one or some nonlinear programming (NLP) solvers. General-purpose optimal control software (e.g. GPOPS-II [?], Muscod-II [?], [?], Acado [?]]) address ii) to iv) but they need to be interfaced with a musculoskeletal modeling module and they do not provide any built-in biomechanics features (physiological cost functions, kinematic constraints, etc.). SCONE is definitely oriented towards biological motion, but it relies on a single shooting transcription, which as shown to be limiting on complex problems. The aforementioned OpenSim Moco, is a welcome initiative that draws its strength from its integration with the widely used OpenSim. However, it faces the following limitations: i) it is interfaced with only one solver (IPOPT), ii) it only implements direct collocation as transcription method, preventing the possibility to use multiple shooting and iii) it is not as flexible as required by the community, since it requires C++ knowledge for the user to develop new features, such as new objective functions.

The objective of the present paper is to introduce Bioptim1 , an open source optimal control software dedicated to musculoskeletal biomechanics. Bioptim is based on C++ code for computational efficiency but the user interface is written in Python for flexibility and ease-of-use. The OCP transcription can be either direct multiple shooting to preserve the possibility of using arbitrarily accurate ODE solvers for the integration, or direct collocation for the user to choose the most efficient formulation to solve its problem. In both cases, the forward dynamics function is fully parallelized over the shooting intervals for more efficiency. Bioptim's core is fully written in CasADi symbolics to benefit from algorithmic differentiation and to exploit CasADi's interface with several non-linear solvers (Ipopt, SNOPT). Moreover, Bioptim is interfaced with the cutting-edge solver Acados, a recent NLP solver dedicated to direct multiple shooting, intended for realtime applications. The purpose of Bioptim is to allow fast and flexible musculoskeletal OCP formulation and solving by providing a framework with a lot of typical biomechanics problems already implemented and customizable.

The paper is organized as follows: first, the design and implementation of Bioptim are described. Next, the versatility and performances of Bioptim are shown through a variety of examples available online.

II. IMPLEMENTATION AND DESIGN

A. Implementation and dependencies

Bioptim is the top layer of a series of libraries (Biorbd: dynamics and musculoskeletal modeling; CasADi: automatic differentiation; Ipopt/Acados: optimization; Bioviz: visualization). Within this software suite, Bioptim's main role is to shape the problem to allow its dependencies to communicate efficiently, while providing an intuitive and flexible interface to the user (Fig. 1). Therefore, it was written in Python for its flexibility and its widespread use among researchers. However, all intensive calculations behind the interface are performed in C/C++, keeping Bioptim both fast and easy to customize.

B. Design

Bioptim shapes and solves optimal control problems whose two required entries are a model (from Biorbd) and an OCP. The model (editable text file) contains the geometrical characteristics and the segment inertial parameters as well as optional elements, namely, the markers, the actuators of the model (muscles and joint torques possibly with torque/angle/velocity relationships) as well as bounds on joint kinematics and torques. It also allows the user to design or import meshes for visualization purposes. The OCP consists in a combination of nonlinear problems (NLPs) that allows for the formulation of multi-staged OCPs. Each NLP has the following programming attributes: 1) a Dynamics, 2) a set of ObjectiveFunctions, 3) a set of Constraint, 4) Bounds, 5) a number of discretization intervals and the duration of the problem and 6) a set of InitialGuess. Based on these inputs, Bioptim properly sets up the transcription of the OCP, using either direct multiple shooting with several integrators (RK4, RK8, implicit RK and those available in CVODES [?]) or Legendre/Radau direct collocation up to 9 th order with dynamics satisfied at each collocation point. Then the problem is shaped up to feed the chosen nonlinear solver (Ipopt or Acados). Next, we develop the different attributes of each NLP:

1) Dynamics: The Dynamics defines which variables are states (x), controls (u) and parameters (p), the latter being time-independent. Essentially, it implements the ODE governing the system dynamics:

ẋ = f (x, u, p).
(1)

Several Dynamics are already implemented in Bioptim, among which the controls (piece-wise constant/linear, linear continuous) can be muscle excitations, muscle activations and/or joint torques. The states can be muscle activations and/or joint kinematics. The dynamics can include soft (Hunt-Crossley model, [?]) or rigid contact points, external forces, etc. Even if these dynamics types exhaustively span the current usages in biomechanics, a custom dynamics type is also preimplemented to easily customize problems. Bioptim offers the possibility to scale the parameters, if their optimal values are known to differ by several orders of magnitude, which could jeopardize the convergence of the problem.

2) Objective function set: In line with the optimal control formalism, there are two main types of ObjectiveFunctions, namely Lagrange and Mayer. Lagrange types are running objectives, integrated over the NLP duration. Mayer types are time-specific objectives. Classically, they correspond to a terminal objective, but to be more versatile, they can be defined at any instant in Bioptim.

ObjectiveFunctions can depend on any of the optimization variables, i.e. the controls, the states, the parameters and the duration of the problem. Plenty objective function types are already implemented in Bioptim, among which tracking / minimizing, on states / controls / markers / contact forces / problem duration, etc. Should one go missing, a custom objective type is also possible to define.

When declaring the desired list of objective functions for a given NLP, each objective function type is associated with a weight, and the user can choose on which components of the vector variables the objective must apply. If applicable (for tracking objective functions mainly), the user must also specify the numerical target of the objective.

3) Constraint set: Classically, constraints are hard penalties of the optimization problem, i.e., a solution will not be considered optimal, unless all constraints (equality or inequality) are met. The Constraint class contains a variety of implemented constraints. Some of them are specific functions, commonly useful in biomechanical problems (e.g. nonslipping contact point, non-linear bounds on torque depending on the state, etc.), the others have their equivalent in the ObjectiveFunction class. Should one go missing, a custom constraint type is also possible to define.

4) Bounds: Essentially, the Bounds are constraints directly related to the states, the controls and the parameters. They are useful to define model-related constraints such as kinematic, torque or muscle excitation / activation limits.

5) Shooting points and problem duration: In direct approaches, the total duration of the problem is divided into smaller intervals whose initial values are called shooting or collocation points. In Bioptim, the user is asked to define a number of intervals and a problem duration, per phase. Possibly, the problem duration can be part of the optimization variables, allowing for, e.g., minimal time formulations.

6) Initial guesses:

The user can provide an InitialGuess for all the optimization variables, at the beginning of each interval. This feature aims at providing prior information to the solver. Several InterpolationTypes are implemented (constant, linear, spline, each point, etc.), to quickly let the user define the initial guesses. A custom InterpolationType is also possible to implement.

C. Graphical digest

A graphical digest of the OCP can be generated in Bioptim. It consists of a .pdf file auto-generated with Graphviz [?]. It displays the structure of the created OCP as well as all the information for the user to check at a glance the formulation of the problem (Fig. 11).

-⋆ - - 1e -5 1e -11 1e -7 1e -3 1e -3 1e -3 1e -3 - 1e -3 1e -3 S-S rot. error ( • ) or Time to S-S rot. error ≥ 10 • (s) + 0.08 • 3.6 • 0.1 • 1.16 s 1e -5 • 1e -1 • 5
.64 s 5.05 s 0.12 s 0.14 s 1.14 s 0.14 s 0.14 s

Inside an example, green, orange and red cells compare the performances of the approaches. Grey stands for not comparable.

† In DMS the integrator is RK4. The DC is based on 4th order Legendre polynomial. * Currently, not all types of OCP can be solved using Acados. ‡ Tests conducted on a AMD Ryzen 9 CPU @ 3.40GHz × 8 ⋆ -stands for non applicable. + The single shooting (S-S) state trajectory is obtained by forwardly integrating the initial state with the optimized control inputs using Scipy's RK45 integrator.

The S-S error was computed as the average error between the optimized state vector and the single shooting one at each time.

Due to the precision of the integrator used for computing S-S trajectories and to the length of the kinematic chains, the rotational part of the error tends to diverge.

For errors ≥ 10 • , instead of reporting the error in • at the end of the motion, the time at which S-S rotation error becomes ≥ 10 • is reported.

The translational part of this error is reported mm.

III. EXAMPLES

In Appendix, six examples are presented to illustrate the versatility of Bioptim and give a practical overview on how to use its main features. They include an excitation-driven pointing task and a moving horizon estimation of muscle forces with arm models, twisting somersault, walking gait and jumping with whole-body models and a stabilization task with a pendulum-on-a-spring model. The performances (number of solver epochs, convergence time, integration error, optimized objective, influence of the transcription method) of each of these examples are summarized in Tab. I. When possible, problems were solved with both Ipopt and Acados. The source code of these example is available here. Beyond what is presented in this paper, every feature of Bioptim is thoroughly illustrated by additional examples available online, in the getting started folder of the project (parameter optimization, custom objects, multi-phase problems, etc.).

IV. DISCUSSION

The purpose of Bioptim is to solve a variety of biomechanical OCPs with minimal user effort and high performances in terms of computational time. The main features illustrated by the six provided examples are (Tab. I):

• the possibility to use torque-or muscle-driven models (and their combinations); • a variety of ready-to-use cost functions, constraints and dynamics (with and without contacts)... • ... easily customizable in Python when required by the user; • the possibility to solve advanced OCPs (possibly multiphase) in a few seconds or minutes, that previously took us hours; • the possibility to use, in a same software direct collocation or direct multiple shooting; • the interface with two different NLP solvers.

A. Choose your transcription method!

As illustrated by the lack of consensus emerging from Tab.I, the debate remains about the performances of direct collocations versus direct multiple shooting [?], [?]. As a consequence, Bioptim is a powerful tool to investigate this question. Concerning the integration method in multiple shooting, either internally or via Acados, several schemes are implemented (RK4, RK8, implicit RK). While IRK coupled with Acados showed better convergence in our experience with hard problems, RK4 showed to be a good speed/robustness tradeoff in most of the cases. In our experience with Bioptim, direct multiple shooting is not a limitation to the performances (cost value and time to convergence). The variety of transcription possibilities (methods, integration schemes), will allow us, in future studies, to provide useful recommendations for solving biomechanical OCPs and to complement the findings of [?].

B. Automatic differentiation

One of the reasons explaining the performances of Bioptim is the rewriting of the core software, RBDL [?] and Biorbd implementing the dynamics, into CasADi symbolics, as in [?], to automatically provide the exact Jacobians and Hessians of the resulting NLP. The gain in accuracy for the calculation of derivatives leads to shorter convergence times (due to much less iterations) and to optimal solutions reached with lower tolerances. This last aspect must be emphasized for complex motions (fast, highly dynamics ones), because, for instance when using Ipopt, an optimal solution obtained with an unscaled convergence criterion of 10 -2 is very unlikely to be dynamically sound; In other words, it would diverge when forwardly integrating the controls in a single-shooting manner. A lower tolerance (10 -6 or 10 -8 ), which is only reachable with exact derivatives-for most of OCPs in biomechanics-, is expected to lead to better forward dynamics results.

C. Python based, but fast!

Bioptim was thought as an interface, and was therefore written in Python to allow the user to easily combine existing cost functions or constraints and self-implemented ones, to switch from one solver to another, etc. We believe this feature to be of importance given that the biomechanics community is mainly composed of software users rather than developers. Therefore, providing a custom interface in Python rather than in C++, was a driving objective of our work to facilitate a rapid appropriation by the community. Since flexibility and ease-of-use should not compromise the performances, the integration is multi-threaded (openMP-based) and all the inside computations are expressed as C++ CasADi graphs, interfaced with C++ NLP solvers. These graphs can either be built in casadi.MX() or casadi.SX(). The latter requires more RAM for building the problem but is faster to solve. While both may be used with Ipopt, Acados is only compatible with casadi.SX(). By leveraging the speed of casadi.SX() graphs, we were able to estimate muscle forces in real time using Acados on a standard laptop (Appendix E). For a more in-depth analysis of the real-time estimation capabilities of Bioptim, see [?]. Alongside with the 3D visualizer Bioviz that animates the solution, Bioptim proposes a series of onlinegenerated figures, inspired by the real-time graphics from Muscod-II [?], [?], to visualize the optimized variables at each iteration of the solver. This is made with minimal computational cost thanks to the multiprocessing Python toolbox. Our implementation leverages the Python pickle library for easily saving and loading OCPs for, e.g., post-processing analysis. Finally, every layer (integration, optimization, visualization) of Bioptim is optimized to be flexible and fast.

D. Fast vs robust NLP solvers

Fast solvers, such as Acados, offer the opportunity to use multi-start approaches on complex problems, to circumvent the obstacle of local minima [?], [?]. It also allows to get meaningful initial solutions from simpler problems, for guiding the resolution of the harder problems. On the other hand, robust solvers, such as Ipopt, are convenient when the user lacks information about the sought solutions and thus cannot guide the solver through a good initial guess. For biomechanics applications, the complementary characteristics of the interfaced solvers is a really useful tool. Moreover, Bioptim's full compatibility with CasADi provides the opportunity to use any solver already interfaced with it, including third-party software such as SNOPT, WORHP [?] and KNITRO [?] (not tested yet).

E. Multiphase

Biomechanics studies often face changing dynamics or objective functions due to the loss or gain of contacts or timevarying biomechanical tasks. When tracking such a motion or trying to predict it, these changes translate into multiphase OCP. This is one of the reported drawbacks of OpenSim Moco, which does not provide this feature yet. Bioptim, however, is able to handle multiphase OCPs, although they can currently only be solved with Ipopt (Appendix D and F).

F. From constraints to objectives: easy problem relaxation

As stated in Sec. II.B, there exists a correspondence between most of the pre-implemented Constraints and ObjectiveFunctions. This is intended to allow for easy relaxation when the problem is reluctant to converge. For instance, when a biomechanical task requires the final configuration of the model to be enforced (reaching, cyclic motions, sports, etc.), one should first use a Constraint (e.g., TRACK STATE). If the convergence is challenging, just turning this constraint into its namesake Mayer ObjectiveFunction, with a heavy weight, should help the solver.

G. Limitations

Bioptim is already a mature solution for solving biomechanical OCP. However some limitations should be raised. First, it is based on Biorbd which is not as advanced as OpenSim or AnyBody (AnyBody Technology) in terms of biomechanical features and audience. Nevertheless, Biorbd is actively maintained, fast and CasADi-compatible for automatic differentiation. The variety of proposed examples highlighted simple to advanced models. Even if defining a new model was made straightforward thanks to the .bioMod file format, biorbd does not include a GUI for building models. Some Opensim models can be translated into .bioMod but Biorbd curently only supports one single wrapping object per muscle and rigid tendons. As seen in [?], wrapping objects are rare in biomechanical OCPs due to the computational cost and required optimization when a line of action is in contact with more than one object, which compromises automatic differentiation. Via-points and pre-processed moment arms [?] (to be expressed as polynomial functions of crossed DoFs) are often preferred.

H. Future directions

Bioptim v2.1.1, was released in November 2021, with all the features presented in this communication. Some improvements are expected in a near future. First, a graphical model builder is planned in Biorbd, to easily generate .bioMod files. Then, soft contacts are implemented, but they still yield to some instabilities that must be worked on. Also, fatigue models are to be included in Bioptim, to predict adapted motor strategies for long or demanding motions. The implementation of muscle-tendon equilibrium is planned for fast movements or those with large ranges of motions. It will require an additional optimization step to achieve the equilibrium as done in CEINMS [?] or the addition of muscle lengths as state variables, as in [?]. Moreover, an effort will be made to extend the compatibility of Acados with all the features of Bioptim (multiphase, nonlinear constraints, etc.). We also plan to add an inverse optimal control module to Bioptim and muscle synergy dynamics to improve motion predictions [?]. Finally, Bioptim needs to be properly benchmarked against other existing solutions to demonstrate its accuracy and rapidness.

ACKNOWLEDGMENT

Bioptim acts as a catalyst in our group and several students contributed to this library. Thank you to A. Dang, T. Gousselot, P. Puchaud, V. Thiron, L. Vayssac, A. Venne and P. Wegiel.

APPENDIX

In this appendix, six examples are presented to illustrate the versatility of Bioptim and give a practical overview on how to use its main features. The settings and performances (convergence time, single shooting integration error, optimized objective) of each OCP are summarized in Tab. I. When possible, problems were solved with both Ipopt and Acados. In the following, bold symbols denote vectors and hat ones ( * ) denote reference or tracked quantities. In these examples, when muscles are modeled, they are of Thelen's type [?].

A. Muscle excitation driven pointing task

In this first example, the goal was to achieve a muscle excitation driven pointing task using a 2-DoF arm model with six muscle elements. In addition to muscle-induced torques, pure joint torques were added to compensate for the model weaknesses. The main term (highest weight) of the objective function (Eq. 2) is a Mayer objective, corresponding to the pointing tasks at the final node, to superimpose two markers, the first one (m ) fixed in the ulna system of coordinates and the second one (m * s ) fixed in the scene. The three Lagrange terms were added for control (muscle excitations e and joint torques τ ) and state (x) regularization:

J = ω 1 ∥m u (T ) -m * s ∥ 2 TRACK MARKERS T t=0 ∥e∥ 2 MIN EXCITATION + ∥τ ∥ 2 MIN TORQUE + ∥x∥ 2 MIN STATE dt, (2)
where T = 2 s is the duration of the motion, and ω 1 = 1e5. The problem was discretized using 200 shooting nodes with a 5-steps Runge-Kutta-4 (RK4) integration in-between. The problem was solved using Ipopt (with exact Hessian computations) and Acados (with a Gauss-Newton approximation of the Hessian) resulting in two very close solutions. Acados was about 50 times faster than Ipopt and was better at enforcing the continuity constraints (as shown by the single shooting error in Tab. I). Ipopt however ended up with a smaller optimized objective (20.8 vs 23.2), leading to a more optimal solution than Acados. Superimposed snapshots of the optimal motion found with Acados are displayed in Fig. 3. It is worth mentioning that for the purpose of this illustration, no constraint was given on the shoulder range of motion to ensure physiological muscle trajectories.

B. Quaternion-based twisting somersault

In this example of acrobatic sports biomechanics, the goal was to maximize the twist rotation (ϕ) of an 8-DoF model in a backward somersault. It illustrates Bioptim's ability to handle quaternionic representations of rotations.

The model was composed of a 6-DoF root segment and two 1-DoF torque actuated shoulder joints. Two different numerical descriptions of the root segment rotations were used: Euler angles and quaternions. The objective function was as follows:

J = T 0 ω 1 φ MIN TWIST + ω 2 ∥τ ∥ 2 MIN TORQUE dt,
(3) with ω 1 = -1 (resulting in the maximization of the first term) and ω 2 = 10 -6 , T is the duration of the movement and τ is the torque control vector. The first term of the objective function (Eq. 3) corresponds to maximizing the change in twist rotation and the second term is for control regularization. The movement lasted for approximately 1 second and was discretized with 100 shooting nodes, a kinogram is presented in Fig. 5. The optimal kinematics were different for the two types of models (Fig. 4) because of the presence of local minima. However, both models take profits of a common biomechanical strategy (i.e. tilting the body to bring closer together the twist axis and the angular momentum vector), leading to a close performance (2.5 twists), which highlights the equivalence of the two representations. Euler angles have the advantage to be easily interpretable, but they suffer from the loss of a DoF at the gimbal lock (leading to numerical instabilities). The use of a quaternion-based representation tackles this numerical stability issue when a joint is free to rotate on a three-dimensional range of motion. Quaternion's integration must be handled with care [?]. Indeed, when representing orientations, quaternions must be unitary and thus belong to a constrained manifold (namely, the unit 3-sphere S 3 ). However, classical numerical integration schemes such as Runge-Kutta methods treat unit quaternions as if they were arbitrarily defined in R 4 . To overcome this challenge, Bioptim performs a normalization after each Runge-Kutta iteration to project non-unitary quaternions onto S 3 .

C. Pendulum on a spring

This example is presented to introduce Bioptim's ability to use external forces. The goal was to hold the position of a 1 kg mass hanging on a linear spring attached to the ground. A 0.2 m-long pendulum weighting 10 kg was attached to the mass and free to rotate in one dimension (Fig. 6). In addition to the spring force, the mass was actuated by a vertical external force, τ , (e.g., something pulling on it) while the pendulum rotation was passive. The system therefore comprised two DoFs, the mass position (q m ) and the pendulum angle (q p ) and one control input, the vertical external force pulling on the mass (τ ). The spring force F s was:

F s = -k * q m , (4)
with k the spring stiffness constant. The OCP was composed of two phases each lasting for 5 s, with 50 shooting nodes. In the first phase, no objective function was minimized and τ was constrained to be 0, letting the mass oscillating freely. Then, in the second phase, a cost function (Eq. 5) was minimized, to enforce a reference position q * m of the mass. This objective function, exclusively composed of Lagrange terms, was formulated as follows:

J = T T /2 (q m -q * m ) 2 TRACK STATE + ω 1 τ 2 dt MIN TORQUE dt, (5)
with q * m = -0.5 m and ω 1 = 10 -6 and T is the duration of the movement. The first term of the objective function (Eq. 5) acts as a position controller for the mass. The second was added for control regularization.

During the first phase, the mass is passively oscillating around its stationary position due to the spring force (Fig. 7). At the beginning of the second phase, when an additional external force acts on the mass, it stabilizes around the targeted position. The standard deviation between the position and the targeted position is 0.04 m. This example highlights the possibility of using optimal control to find activation patterns compensating for external passive forces (e.g., ortheses flexibility, contact surface deformation, interaction between two models, etc.).

D. Multiphase activation driven walking cycle

This example is presented to introduce Bioptim's ability to deal with a multiphase locomotion estimation problem including muscle actuation and contact forces. The goal was to estimate muscle activations by tracking markers trajectories and ground reaction forces. The model was a 3D leg with 13 DoFs (6-DoF pelvis, 3-DoF hip, 1-DoF knee, 2-DoF ankle and 1-DoF toes), driven by 19 muscle activations and residual joint torques to compensate for potential muscle actuation weaknesses. The gait cycle was defined from the first heel strike to the end of the swing phase discretized into 90 shooting intervals. To approximate the natural rolling of the foot, the stance was divided into three phases (heel, flatfoot and forefoot contacts) of fixed duration deduced from experimental force platform data and markers position (0.05, 0.36 and 0.16 s). The swing phase lasted 0.38 s. The interaction between the ground and the foot was modeled using a fourcontact points model located at the heel and the forefoot (first, fifth metatarsi and hallux). The optimization problem consisted in minimizing the errors between predicted (m) and reference (m * ) markers trajectories, predicted (F , M) and reference (F * , M * ) ground reaction forces at all contact points. A regularization term on muscle activations (a) was also added (least-activations) as well as a penalization term on the residual torques (τ ):

J = 4 i=1   Ti Ti-1 ω1(∥m -m * ∥ 2 ) TRACK MARKERS + αω2(∥F -F * ∥ 2 ) TRACK FORCES + ω3∥a∥ 2 MIN ACTIVATION + ∥τ ∥ 2 MIN TORQUE dt , (6)
where ω 1 = 10 4 , ω 2 = 10 -2 , ω 3 = 10 are weighting factors, T 0 = 0 and T i , i ∈ [1, 2, 3, 4], are the final time of the i th phase. Ground reaction forces and moments were only tracked during the stance phase, hence α = 0 during the swing phase and α = 1 otherwise.

Non-slipping (NON SLIPPING) and unilateral contact force (CONTACT FORCE) constraints were added to prevent the foot from slipping and pulling from the ground. In between phases, the use of the PhaseTransition.IMPACT state transition allowed to represent the gain or loss of contact(s) in the dynamics (e.g., [?] swing phase to heel strike ).

Tracking experimental data allowed to reproduce leg motion during the walking cycle (Fig. 8). The root mean square tracking error on markers trajectories was 17 mm (mean errors on pelvic and foot markers were 10.4 mm and 10.2 mm, respectively). Concerning ground reaction forces tracking, the root mean square error was 7 N. During the stance phase, Gluteal muscles were mainly activated during the stance phase especially during initial contact (1 %), followed by the Vastus Medialis for the loading response (10 %). We found light hamstrings activity during the early stance and terminal swing.

The transition from stance to swing (60 % -70 %) was highly actuated by hip flexors (Iliopsoas and Rectus Femoris) and leg muscles (Gastocnemius Medialis). Activation of the Tibialis anterior occurred at the beginning and end of the cycle, probably to adjust foot position (Fig. 8). These results were similar to the characteristic average activity patterns of the lower limb muscles during locomotion described in [?].

E. Moving Horizon Estimation of Shoulder Elevation

This example is presented to introduce Bioptim's ability to provide real-time estimation of biomechanical variables. The goal was to perform a real-time estimation of dynamically consistent joint kinematics and muscle forces, using a moving horizon estimation (MHE) approach (i.e. an optimization approach that uses a series of measurements observed over time). A shoulder elevation motion was performed with a 4-DoF (q) arm actuated by 19 Hill-type muscle elements. The control inputs of the model were the muscle activations (a). The MHE implementation consists in splitting the OCP into a succession of smaller ones for processing fixed-size subsets of the tracking data moving forward in time. Each time one subproblem is solved, a new measurement is added, the oldest one is discarded and a new subproblem is defined. Due to their similarities, the solution of the previous OCP is a good initial guess to the new one. The dynamical consistency of the final solution is enforced by continuity constraints on the initial state. Each objective function (Eq. 7) was written as the sum of three terms: tracking reference joint angles (q * ), states and muscle activations regularizations (i.e., least-square criteria):

J = t+tmhe t ω 1 (∥q -q * ∥ 2 ) TRACK STATE + ω 2 ∥q∥ 2 MIN STATE + ω 3 ∥a∥ 2 MIN ACTIVATION dt, (7)
where ω 1 = 10 3 , ω 2 = 10, ω 3 = 10 2 and t mhe is duration of each sub-problem.

In this example, reference data of an 8 s series of four arm elevations were generated at 100 Hz, by computer simulation. A centered Gaussian noise (mean = 0, std = 0.04) was added to q * , to simulate experimental-like joints angle measurements. Using a window size of 7 nodes (i.e., 210 ms), the estimator ran at about 33 Hz (one in three reference data frame was sent to the estimator to simulate experimental-like conditions), i.e., two and half times faster than standard biofeedback (13 Hz, [?]). The MHE was able to forecast the movement kinematics with a root mean square error of 1.3 ± 0.7 • while providing a realistic estimation of muscle forces close to the ground truth with a root mean square error of 11.1 ± 14.9 N (Fig. 9).

F. Multiphase vertical jumper

This example was designed to introduce Bioptim's ability to reduce the number of degrees-of-freedom (DoF) of a model via the BiMapping feature, to account for nonlinear boundaries on the controls, and to solve complex multiphase OCP. Two phases were used to describe the dynamics of the push-off phase of the jump: flat foot (two floor contacts) and then toe only (one contact) 2 . A pseudo-2D full-body symmetric model consisting of 3 DoFs at the pelvis (forward and upward translations, tranverse rotation), 1 DoF at the upper limb (shoulder flexion), and 3 DoFs at the lower limb (hip, knee and ankle flexion) was used. Since this is a full-body model, the root segment (i.e., the pelvis) was left uncontrolled, reducing the number of control variables to four, namely the shoulder, hip, knee and ankle flexions. The objective function with the most important weight was a Mayer objective computed at the end of the push-off phase consisting in maximizing the jump height (h) from the free fall equations applied to the center of mass.

J = ω h h MIN PREDICTED COM HEIGHT + 5 i=1 ω t (T i -T i-1 ) MIN TIME dt (8)
where T i with i ∈ [1, 2] are the final times of the i th phase respectively, and T 0 = 0; ω h = -100 is the weight of the jump height term defined negative to maximize it and ω t = 0.1 is the weight put on time minimization; Joint angles were bounded to human-like limits. The first node of the first phase was constrained such that the avatar 2 A contact is defined as a point where forces are applied to cancel its acceleration. stood with its knee slightly flexed and its arms horizontal. Joint velocities were arbitrarily bounded to human-like limits. Joint torques were bounded with nonlinear torque/angle/velocity relashionships measured on a high-level athlete using an isokinetic dynamometer [?]. Non slipping (NON SLIPPING) and unilateral (CONTACT FORCE) contact force constraints were added to prevent the contact points from slipping and pulling on the ground. During the second phase, the heels had to remain over the floor. To speed-up the convergence with Ipopt, the problem was first solved using a BFGS Hessian approximation. Then, starting from this first solution, the problem was re-optimized, with the exact-Hessian computations.

The optimized solution was obtained in 75 + 33 (resp. 175 + 26) iterations for the direct multiple shooting (resp. direct collocation), resulting in a 1.98 m jump height. The solution reproduced a human proximo-distal strategy (Fig 10 ), i.e., activating large segments first (for instance the torso) and sequentially adding more distal segments, consequently ending up with the feet.

G. Graphical Digest

In the following, an example of auto-generated graphical digest is presented. It sums up the OCP of Ex.C, highlighting the content of each phase of the problem (including, e.g., the integration method),and the way the transition between them is handled.

Fig. 1 .

 1 Fig. 1. Bioptim dependencies flowchart. The red-boxed software are developed by the S2M team. The Bioptim part is further detailed in Fig. 2.

Fig. 2 .

 2 Fig. 2. Bioptim design flowchart. Red boxes correspond to objects that must be filled in by the user. Red-dashed boxes correspond to preimplemented objects already available to the users. ⊗ stands for easily customizable objects.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Snapshots of an optimized activation-driven pointing task with Acados. The arm starts facing upwards in left hand part of the picture and ends facing downwards in the right hand part. The marker fixed on the ulna head is depicted in blue and the scene-fixed target marker is depicted in red. Red lines show the lines of actions of the muscles.

Fig. 5 .

 5 Fig. 5. Snapshots of maximally twisting somersaults driven by shoulder torque actuators and a free base whose rotation is expressed with Euler angles.

Fig. 7 .

 7 Fig. 7. Two-phases kinematics of the mass-pendulum-spring system. Gray dashed lines show the phase transition, blue lines are related to the mass (position velocity and external force acting on it), red lines are related to the pendulum (position and velocity) and the green line depicts the spring force.

Fig. 8 .

 8 Fig. 8. Snapshots of a walking gait cycle driven by muscles activation with histogram of muscle activations below. On top, the percentage of the gait cycle are displayed. The red lines represent muscles lines of action and the blue points depict the tracked markers. The activation of the Gluteus Maximus is the mean of its three parts and the Hamstring is the mean activation of the Semimembranous, Semitendinous and Biceps Femoris.

Fig. 9 .

 9 Fig. 9. Ex. E. Top row -Real-time estimated joint angles (blue), ground-truth joint angles (orange) and tracked noisy joint angles (green) for a cyclic motion of the arm. Bottom row -Real-time estimated muscle forces (blue) and ground truth muscle forces (orange) for the same motion. Only four muscles with significative action (peak forces > 15 N), on the two selected DoFs, are shown. Muscle abbreviations stand for (from left to right): Triceps Long head, Deltoid Middle, Infraspinatus, Biceps Brachial Short head.

Fig. 10 .

 10 Fig. 10. Snapshots of the push-off phase of a vertical jump (Ex. F). The avatar reproduces a human-like jump movement. The first three positions represent the first phase of the optimization (i.e., heel and toe in contact with the floor) and the fourth position depicts the end of the second phase (i.e., only the toes in contact with the floor)

Fig. 11 .

 11 Fig. 11. Graphical digest of problem C. Benjmain Michaud received the M.Sc. degree in biomechanics in 2012, and has since been working in this field. He is the lead developer of ezc3d and of the biorbd suite. In 2019, he has received the prestigious Vanier Scholarship to begin its PhD on predictive method to reduce muscle fatigue of violinists within the Laboratory of Simulation and Movement Modeling at Universit é de Montr éal, Canada.

TABLE I OVERVIEW

 I OF THE COMPUTATIONAL RESULTS AND DYNAMIC CONSISTENCY FOR SIX EXAMPLES THAT INCLUDE A VARIETY OF FEATURES OF Bioptim.

			A. Pointing			B. Somersault		C. Pendulum	D. Walking	E. MHE	F. Jumping
					Euler	Quat.
	Transcription †		DMS	DC	DMS DC	DMS	DMS DC	DMS DC	DMS	DMS	DC
	Solver *	Ipopt	Acados	Ipopt	Ipopt	Ipopt	Ipopt	Ipopt	Ipopt	Ipopt	Ipopt	Acados	Ipopt	Ipopt
	# NLP iterations	39	58	40	217	41	186	116 129	832 352	-	75 + 33 175 + 26
	Optimized cost	21.6	277	21.5	-21.5	-6.3	-21.8	5e -3 7e -3	94.5	95.4	-	-198.3	-198.3
	Time to convergence ‡	27 s 2.7 s	5.9 s	112 s 20.7 s	90 s	25.8 s 27.5 s	8.6 h 1.7 h	6.2 s	163 s 76 s
	S-S trans. error (mm) +

link -DOI: 10.5281/zenodo.5711443

This work was partly funded by a scholarship of the Vanier program (BM), the Canada First Research Excellence Fund via the Trans-MedTech Institute (FB) and the NSERC Discovey Programme (MB).

* B. Michaud and F. Bailly have contributed equally to this work and share first authorship.


Abstract : A chemo-and stereoselective palladium-catalyzed amination of silylated butenediol dicarbonates has allowed for the introduction of a glycine moiety to obtain a desired functionalized epoxysilane. A stereoselective aldolization then delivered the piperidine ring which may be used as a precursor for the synthesis of a variety of polyhydroxylated azasugars. This efficient approach has been illustrated by the synthesis of 1-deoxymannojirimycin including a stereoselective reduction with LAH and a Tamao-Fleming oxidation of a C-SiMe2Ph bond.

For nearly half a century, glycosidase inhibitors have been the subject of intense interest. 1 Polyhydroxylated derivatives of nitrogen heterocycles are one class of active compounds that can be used as potential drugs to treat diabetes, hepatitis, or various cancers. 2 At physiological pH, their ammonium salts have been recognized to act as transition-state analogues of the enzymes. Moreover, as has been pointed out by Bols, 3 "The best glycosidase inhibitors are compounds that have a nitrogen in place of the exo-or endocyclic oxygen or anomeric carbon...It is more important to mimic the charge development in the transition state than the shape." Some polyhydroxylated piperidines which have been well studied for their inhibiting properties are 1-deoxymannojirimycin (1), 1deoxynojirimycin ( 2), and fagomine (3) (Figure 1). Substituted piperidines, in general, 4 and polyhydroxylated derivatives, 5 more precisely, have been the targets of a large number of synthetic approaches 6 due to their potential as therapeutic agents and the need to find efficient access to biologically active analogues. Herein, we propose a general approach for the synthesis of 1-deoxyazasugars as well as some silylated derivatives. 7 Scheme 1 illustrates our retrosynthetic analysis. Access to the desired polyhydroxylated piperidines could be achieved in a divergent way from the bicyclic system 4. The nitrogen heterocycle bearing the epoxysilane function was considered to be of great interest, as it could be functionalized following selected reactions such as nucleophilic attacks on the epoxysilane 8 or Brook-type rearrangements, 9 with or without migration of the silicon group. 10 For example, from 4, chemoselective reduction of the ester and epoxysilane functionalities should deliver the silyl analogue of DMJ. Further reaction utilizing a Tamao-Fleming oxidation of the C-Si bond should yield the protected form of 1. 11 On the other hand, direct oxidation of the epoxysilane 12 should give a keto alcohol, which can be stereoselectively reduced to yield DNJ 2. Our synthetic approach to the piperidine skeleton would be based on an intramolecular aldol condensation of the acyclic precursor 5. [START_REF] Danieli | For related intramolecular aldolisations[END_REF] Previously, we reported a strategy for the enantioselective preparation of silylated epoxy cyclopentanols via a stereoselective aldolization to build the five-membered cycle. 14 We expected the same high level of stereoselectivity in 5 due, in part, to the presence of the three-membered ring.

Scheme 1. Retrosynthetic Analysis of Polyhydroxypiperidines

Finally, the preparation of the desired allylamine 6 from the easily accessible silylated bis-allylic derivatives 7 would nicely illustrated the chemo-and stereoselective silicon directed palladiumcatalyzed alkylation we have recently published. 15 Our synthesis began with the preparation of the model diacetate 7a bearing the triethylsilyl group.

We found it to react with N-tosylglycine methyl ester in the presence of palladium catalyst and i-PrOH as solvent (Scheme 2). [START_REF] Thorimbert | For the reactivity of silylated acetate alcohols with benzylamine[END_REF]17 Unfortunately, due to the low reactivity of the intermediate π-allylic palladium complex, a competing β-elimination also occurred to give the corresponding diene 9a in 40% yield. We therefore prepared a more reactive precursor 7b having two allylic carbonate functionalities. This turned out to be advantageous, the desired product 8b being produced in 75% yield with only a trace amount of the diene 9b. Careful control of the reaction temperature was required for good reproducibility. We noticed similar reactivity with the bis-carbonate 7c, bearing a functionalizable Me2PhSi group. In this instance, we were able to isolate the allylic amine 8c in 71% yield as a 95/5 mixture of two stereoisomers.

Scheme 2. Preparation and Palladium-Catalyzed Amination of Silylated Allylic Derivatives 7

Treatment of 8c with a catalytic amount of potassium carbonate in MeOH allowed for a smooth deprotection to provide the corresponding allylic alcohol 6 in excellent yield (Scheme 3). Epoxidation was next performed in a classic manner with 2 equiv of m-CPBA. The epoxy alcohol 10, so obtained was isolated in 75% yield after purification. However, the crude product was sufficiently pure to be used directly for the next step.

Scheme 3. Preparation of the Piperidine Skeleton

Based on some of our past work, we used the iodoxybenzoic acid (IBX) to oxidize the primary alcohol. We also tested SIBX, 18 a stabilized form of IBX, and obtained the desired aldehyde 5 in a comparable yield of 65%. Significantly, these three consecutive steps could be performed without purification of the intermediates to give 5 from 8c in 85% yield. The acyclic epoxy aldehyde 5 reacted smoothly in the presence of DBU to give the piperidine ring in 92% yield. [START_REF] Keynes | For related intermolecular alkylations or aldolizations of amino esters[END_REF] This aldolization proved highly stereoselective, giving mainly one of the four possible diastereomers (Scheme 3). As expected, the trans relationship between the oxirane and the created hydroxyl was totally controlled. 14 The two diastereomers 4 and 11, obtained in an 85/15 ratio, corresponded to the two epimers at the C-α of the cyclic amino-ester. In the major product 4, the ester functionality was trans relative to the hydroxyl group. Attempts to improve the diastereoselectivity by lowering the temperature, with or without catalytic Yb(OTf)3, were unsuccessful. The two diastereomers 4 and 11 could be separated by careful SiO2 flash chromatography. However, we reasoned that under our basic conditions 11 could give the corresponding dehydroamino ester by dehydration. Indeed, from 5 after a reaction time of 12 h, we obtained an 85/15 mixture of 4 and the expected dehydroamino ester deriving from 11. This allowed for an easy isolation of pure 4 in 78% yield.

With the piperidine ring system in hand, we attempted reduction of both the ester and the epoxide functionalities. At first, the use of an excess of lithium aluminum hydride reduced both the ester and the epoxide in moderate yield (46%). We reasoned that the free hydroxyl might first be chelating to the reducing agent, thus disminishing the regioselectivity of the reduction of the epoxysilyl function. Thus, we protected the alcohol 4 as its TBDMS ether 12 in 96% yield. To our delight, this protection proved to be beneficial for the reduction. After one night in Et2O at rt, we observed the clean formation of two compounds. Direct conversion of the crude mixture to the corresponding acetates greatly simplified the purification (Scheme 4).

Scheme 4. Chemoselective Reduction of the Ester and Epoxide and N-Tosyl Functions

We isolated the expected diacetate 13 in 54% yield, as well as 6% of the triacetate 14 resulting from the removal of the N-tosyl protecting group. The formation of 14 was quite surprising in view of past literature. [START_REF] Greene | Protective Groups in Organic Synthesis[END_REF] Indeed, the N-Ts bond is usually cleaved under more drastic conditions such as DIBAL in refluxing toluene, or by the use of various radical anions, such as sodium naphthalemide. Finally, we took advantage of this result and were able to cleanly convert 12 into the triacetate 14 in 78% yield (Scheme 4).

To explain the clean deprotection of the secondary amine, we propose that the ester of 12 is first reduced into an aluminate, which then assists in the cleavage of the N-Ts bond. The resulting alumino amino-alcoholate could form a stable 5-membered cycle A (Figure 2). This proposition is supported by the fact that we never observed N-Ts bond cleavage for 4. For this compound, the aluminum is complexed by the initially free alcohol early before any reductions and finally gives the stable intermediate B. Interestingly, these polyhydroxylated piperidines are orthogonally protected, which could be helpful for further modifications. Finally, after selective removal of the TBS protecting group, we acetylated the intermediate diol to isolate in a quantitative yield, the known penta-acetylated form of (±)-DMJ 17. Our 1 H NMR spectrum was not in accordance with the one described by O'Doherty, but quite similar to the description of Hutchinson. 23 We fully deprotected 17 with 6 N HCl aqueous solution and obtained the expected product in a quantitative yield. The 1 H and 13 C NMR spectra were in perfect accordance with the descriptions of Han. 24 Finally, the stereochemical assignment of the formed compounds have been confirmed by an X-ray structure analysis of 15 (Figure 3). In conclusion, we have reported a complementary, highly selective, and efficient synthesis of polyhydroxylated piperidines that rivals previous preparations (36%, 7 steps from the acyclic precursor 8c). Our approach allows for the differentiation of all the hydroxy groups, which we believe offer a great potential for the elaboration of more complex racemic molecules. A task which is under active investigation in our laboratory, as well as an asymmetric version of this synthetic approach. 26

Figure 1 .

 1 Figure 1. Structure of 1-deoxymannojirimycin (DMJ) (1), 1-deoxynojirimycin (DNJ) (2), and fagomine (3).

Figure 2 .

 2 Figure 2. Proposed intermediates during the reduction of 4 and 12

  25

Figure 3 .

 3 Figure 3. X-ray Analysis of 15.

A Stereoselective Route toward Polyhydoxylated Piperidines. A Total Synthesis of (±)-Deoxymannojirimycin Cécile

  Boglio, Sebastian Stahlke, Serge Thorimbert,* Max Malacria*

	Institut de chimie moléculaire (FR 2769), Laboratoire de chimie organique (UMR CNRS 7611),
	Université Pierre et Marie Curie, 4 place Jussieu, 75005 Paris, France

Acknowledgment

We gratefully acknowledge financial support from UPMC, CNRS, and Institut Universitaire de France of whom M.M. is a senior member. C.B. thanks the Ministère de l'Education for financial support. S.S. from the University of Göttingen, thanks the Erasmus program. We also wish to thank SIMAFEX and Dr. D. Depernet (dominique.depernet@simafex.fr) for a gift of SIBX.