Cocofest

An Open-Source Python Package for Functional Electrical Stimulation (FES) Optimization in Optimal Control.
Supports predictive musculoskeletal simulation driven by FES, moving time horizon, and model identification.
"Prototype today’s FES to power tomorrow’s rehab."

Table of contents

• About
• Installation
• Features
  • Available FES models
  • Musculoskeletal model driven by FES
  • Moving time horizons
  • Identification
  • Initial value problem
  • Summation truncation
• Other
  • Want to contribute?
  • Contributors
  • Citing
  • Cited in
  • Other related projects
  • Acknowledgements

About

Functional electrical stimulation (FES) is a neurorehabilitation technique that promotes motor recovery after neurological injury. By delivering coordinated electrical pulses to targeted muscles, FES elicits functional movements such as walking, reaching, and grasping. Because responses to stimulation vary across individuals and muscle groups, most FES protocols still rely on empirically tuned parameters. These settings can cause over stimulation, early muscle fatigue on-set, and reduce therapeutic gains.

Advanced control approaches like optimal control-driven FES can improve FES rehabilitation efficiency by personalizing stimulation parameters to a specific task and patient. Therefore, we designed Cocofest (Custom Optimal COntrol for Functional Electrical STimulation) an open-source Python package for optimal control-driven FES. Cocofest relies on bioptim, an optimal control program framework for biomechanics. bioptim uses biorbd a biomechanics library, benefits from powerful algorithmic differentiation provided by CasADi and robust solver like Ipopt.

Important

Cocofest as no clinical clearance and should not be used for rehabilitation purposes.
Don't forget to the repository to show your support and help us grow the community!

Installation

Currently, no anaconda installation is available. The installation must be done from the sources.
Cloning the repository is the first step to be able to use the package.

Dependencies

Cocofest relies on several libraries. So carefully follow these steps to get everything installed to use Cocofest.
First, create a new conda environment

conda create -n YOUR_ENV_NAME python=3.11

Then, activate the environment

conda activate YOUR_ENV_NAME

After, install the dependencies

conda install numpy matplotlib pytest casadi biorbd pyorerun -c conda-forge

Finally, install the bioptim setup.py file located in your cocofest/external/bioptim folder

cd <path_up_to_cocofest_file>/external/bioptim
python setup.py install

You are now ready to use Cocofest!

Features

📊 Available FES models

All models are implemented at the muscle actuator level, making them applicable to a wide range of problems regardless of the specific optimal control problem.

Model Name	Citation	Description / Focus
Veltink1992	Veltink, P. H., Chizeck, H. J., Crago, P. E., & El-Bialy, A. (1992). Nonlinear joint angle control for artificially stimulated muscle. IEEE Transactions on Biomedical Engineering, 39(4), 368–380.	Nonlinear control of joint angles via electrical stimulation.
Riener1996	Riener, R., Quintern, J., & Schmidt, G. (1996). Biomechanical model of the human knee evaluated by neuromuscular stimulation. Journal of Biomechanics, 29(9), 1157–1167.	Biomechanical knee model validated using neuromuscular stimulation.
Ding2003	Ding, J., Wexler, A. S., & Binder-Macleod, S. A. (2003). Mathematical models for fatigue minimization during functional electrical stimulation. Journal of Electromyography and Kinesiology, 13(6), 575–588.	Focus on mathematical models for minimising fatigue.
Ding2007	Ding, J., Chou, L. W., Kesar, T. M., et al. (2007). Mathematical model that predicts the force–intensity and force–frequency relationships after spinal cord injuries. Muscle & Nerve, 36(2), 214–222.	Predicts force–intensity and force–frequency responses post-SCI.
Marion2009	Marion, M. S., Wexler, A. S., Hull, M. L., & Binder-Macleod, S. A. (2009). Predicting the effect of muscle length on fatigue during electrical stimulation. Muscle & Nerve, 40(4), 573–581.	Examines muscle length impact on fatigue under stimulation.
Marion2013	Marion, M. S., Wexler, A. S., & Hull, M. L. (2013). Predicting non-isometric fatigue induced by electrical stimulation pulse trains as a function of pulse duration. Journal of NeuroEngineering and Rehabilitation, 10, 1–16.	Predicts non-isometric fatigue based on pulse duration.
Hmed2018	Hmed, A. B., Bakir, T., Garnier, Y. M., Sakly, A., Lepers, R., & Binczak, S. (2018). An approach to a muscle force model with force-pulse amplitude relationship of human quadriceps muscles. Computers in Biology and Medicine, 101, 218-228.	Models the relationship between pulse amplitude and force.

Note

It is possible to implement more FES models into Cocofest. Adventurous enough to code it by yourself, we are looking forward to read your pull request. Feel free to reach out on discord or submit an issue if you need help.

🦴 Musculoskeletal model driven by FES

In conventional Hill-type muscle model, muscle force ($F_m$) is the product of $a$ the muscle activation, $F_{max}$ the maximal isometric muscle force, $f_l$ the force-length, $f_v$ the force-velocity and $f_{pas}$ the passive force-length relationship:

$$F_m(t) = a(t)\, F_{\max}\, f_l(\tilde{l}_m)\, f_v(\tilde{v}_m) + f_{pas}(\tilde{l}_m)$$
`Cocofest` replaces $a(t)$ × $F_{max}$ by the force obtained using [FES models](#available-fes-models). This approach allows motions driven-FES simulations, meanwhile benefiting from musculoskeletal model properties (e.g., muscle insertion, weight, inertial).

Note

Used force-length ($f_l$), force-velocity ($f_v$) and passive force-length ($f_{pas}$) are those published by De Groot et al., (2016). Those relationships can be activated or not when initializing your OCP. Modification to the following file can be done to have more/different relationships.

💻 A short musculoskeletal FES-driven example

The following example displays a reaching task using the Arm26 model driven by the Ding2007 FES model.

$$\begin{aligned} \min_{x(\cdot),\,u(\cdot)} \quad & \int_{0}^{T} \sum_{i=1}^{n} F_{m,i}(t)\,dt \\[4pt] \text{s.t.:}\quad & q_{\text{arm}}(t) \in [-0.5,\, 3.14], && \forall t \in [0,T],\\\ & q_{\text{forearm}}(t) \in [0,\, 3.14], && \forall t \in [0,T],\\\ & \text{(last node)}\;\; \|p_{\text{hand\_marker}}(T) - p_{\text{target\_marker}}\| \le \varepsilon,\\\ & u(t) = \begin{bmatrix} \mathrm{pw}_{1}(t)\\ \mathrm{pw}_{2}(t)\\ \vdots\\ \mathrm{pw}_{n}(t) \end{bmatrix} \;\; \text{(pulse widths per muscle)}. \end{aligned}$$

Figure 1: Motion performed for the reaching task and associated muscle force production.

Note

Solved in 6.7 second, computer with an AMD Ryzen Threadripper PRO 7965WXs x 48 processor.
Additional information: frequency = 30Hz, n_shooting = 30, step = 0.033s, final time = 1s, integration = Collocation radau method, polynomial_order = 3, solver = IPOPT.

You can find more examples of musculoskeletal model driven by FES in the following file.

⏳ Moving time horizons

For longer time span simulation and apprehend muscle fatigue apparition, Cocofest implements moving time horizons (MHE).

💻 A short MHE hand cycling FES-driven example

$$\begin{aligned} \min_{x(\cdot),\,u(\cdot)} \quad & \int_{0}^{T} \sum_{i=1}^{n} F_{m,i}(t)\,dt \\[4pt] \text{s.t.:}\quad & q_{\text{arm}}(t) \in [0,\, 1.5], && \forall\, t \in [0,T],\\\ & q_{\text{forearm}}(t) \in [0.5,\, 2.5], && \forall\, t \in [0,T],\\\ & q_{\text{pedal}}(t) \in [0,\, 6.28], && \forall\, t \in [0,T],\\\ & \text{(first node)}\;\; \left\|\,\mathrm{center}_{\text{wheel}}(t_{0}) - \begin{bmatrix} 0.35 \\ 0 \end{bmatrix}\right\| \le \varepsilon,\\\ & \text{(first node)}\;\; \left\|\,\dot{\mathrm{center}}_{\text{wheel}}(t_{0}) - 0\right\| \le \varepsilon,\\\ & \text{(last node)}\;\; \left\|\,q_{\text{pedal}}(T) - 6.28\right\| \le \varepsilon,\\\ & u(t) = \begin{bmatrix} \mathrm{pw}_{1}(t)\\ \mathrm{pw}_{2}(t)\\ \vdots\\ \mathrm{pw}_{n}(t) \end{bmatrix} \;\; \text{(pulse widths per muscle)}. \end{aligned}$$

Figure 2: Motion performed for the cycling task, muscle force contribution per section inspired by Quittmann et al. (2025) and muscle force production above 10% of the maximal force.

Note

Solved in 1.02 second, computer with an AMD Ryzen Threadripper PRO 7965WXs x 48 processor.
Additional information: frequency = 30Hz, n_shooting = 60, step = 0.033s, final time = 2s, integration = Collocation radau method, polynomial_order = 3, solver = IPOPT, simultaneous turn per optimization = 2.

🎯 Initial value problem

The initial value problem feature enables forward nonlinear dynamic integration to simulate the model’s behavior from given initial state and controls (i.e., series of pulse trains). This also permits comparison between FES models without using optimal control methods.

For that, the IvpFes class is used to build the problem.

from cocofest import IvpFes, DingModelFrequencyWithFatigue

fes_parameters = {"model": DingModelFrequencyWithFatigue(), "n_stim": 10}
ivp_parameters = {"n_shooting": 20, "final_time": 1}

ivp = IvpFes(fes_parameters, ivp_parameters)

result, time = ivp.integrate()

🔎 Identification

To personalize FES models to simulated or experimental force, Cocofest supports model identification using optimal control.

💻 A short model identification example

$$\begin{aligned} \min_{x(\cdot),\,p(\cdot)} \quad & \int_{0}^{T} \bigl(F(t) - F_{\text{sim/exp}}(t)\bigr)\,dt \\[4pt] \text{s.t.:}\quad & u(t) = \begin{bmatrix} \mathrm{pw}_{1}(t)\\ \mathrm{pw}_{2}(t)\\ \vdots\\ \mathrm{pw}_{n}(t) \end{bmatrix} = \begin{bmatrix} \mathrm{pw}_{\text{exp}1}(t)\\ \mathrm{pw}_{\text{exp}2}(t)\\ \vdots\\ \mathrm{pw}_{\text{exp}n}(t) \end{bmatrix}. \end{aligned}$$

Figure 3: Model identification optimization

Note

Solved in 0.343 second, computer with an AMD Ryzen Threadripper PRO 7965WXs x 48 processor.
Additional information: frequency = 33Hz, n_shooting = 66, step = 0.03s, final time = 2s, integration = Runge-Kutta 4, integration steps = 10, solver = IPOPT.

✂️ Summation truncation

Cocofest also incorporates the recent numerical truncation method to speed up convergence. This method limits the number of past stimulations considered in the dynamics to reduce the dependency on time-varying states.

model = ModelMaker.create_model("ding2007", stim_time=stim_time, sum_stim_truncation=10)

Tip

To determine the value to use for sum_stim_truncation, you can refer to Tiago et al. (2025) or Co et al. (2024).

Other

🙌 Want to contribute?

We are always looking for new contributors to help us improve Cocofest.
Feel free to check our contributing guidelines to get started.

Don't know where to start? Issues tagged with "Good first issues" are a great place to begin!

🤝 Contributors

📝 Citing

Cocofest is not yet published.
Meanwhile, if you use Cocofest, please cite the following zenodo link: 10.5281/zenodo.17068808.

📚 Cited in

Note

If you used Cocofest in your research, please let us know by submitting an issue or a pull request to add your publication to this list.

Other related projects

🙏 Acknowledgements

🌱 Funding

   

Logo and assets design