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.

git clone https://github.com/pyomeca/cocofest.git
cd cocofest

The repository contains an environment.yml file with all of Cocofest's dependencies and a pyproject.toml file for the package itself. The setup requires the following 3 commands from the repository root:

conda env create -f environment.yml
conda activate cocofest
pip install -e .

conda env create reads environment.yml and creates an environment named cocofest with every dependency in one go, and pip install -e . registers Cocofest as a regular importable package, so import cocofest works from anywhere (not just from inside the examples folder), without any relative-pathing issue.

If you would rather pick your own environment name or dependency versions, follow these steps instead.
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 casadi biorbd pyorerun scikit-optimize bioptim==3.4.0 -c conda-forge

Finally, from the root of the cloned repository, install Cocofest itself in editable mode

pip install -e .

You are now ready to use Cocofest!

Every script in the examples folder is directly runnable. For instance, the pulse width optimization script under getting_started can be run with:

python examples/getting_started/optimization/pulse_width_optimization.py

See examples/README.md for an index of every example folder and what each script demonstrates.

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.

Cocofest gives access to every model above through the ModelMaker factory, using the following string keys. Each base key also has a _with_fatigue variant (adds fatigue dynamics) and, for Marion2009/Marion2013, a _modified variant (pulse width instead of frequency as control):

from cocofest import ModelMaker

model = ModelMaker.create_model("ding2007_with_fatigue", stim_time=[0, 0.1, 0.2])

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:

Cocofest replaces $a(t)$ × $F_{max}$ by the force obtained using FES models. This approach allows motions driven-FES simulations, meanwhile benefiting from musculoskeletal model properties (e.g., muscle insertion, weight, inertial).

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

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

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

For longer time span simulation and apprehend muscle fatigue apparition, Cocofest implements moving time horizons (MHE) through the FesMhe (single muscle) and FesMheMsk (musculoskeletal) classes. Each window is solved, then the horizon slides forward by re-using the previous stimulation history to keep fatigue state continuous across windows.

from cocofest import DingModelPulseWidthFrequencyWithFatigue, OcpFes, FesMhe

model = DingModelPulseWidthFrequencyWithFatigue(stim_time=[...], sum_stim_truncation=10)
dynamics_options = OcpFes.declare_dynamics_options(...)  # see full example for every argument

mhe = FesMhe(
    bio_model=model,
    dynamics=dynamics_options,
    cycle_len=cycle_len,
    cycle_duration=cycle_duration,
    n_cycles_simultaneous=n_cycles_simultaneous,
    n_cycles_to_advance=1,
    ...,
)


def update_functions(_mhe, cycle_idx, _sol):
    return cycle_idx < n_cycles  # keep sliding the window until n_cycles is reached


sol = mhe.solve_fes_mhe(update_functions, solver=..., total_cycles=n_cycles, cycle_solutions=...)

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.

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

stim_time = list(range(10))  # 10 stimulations, one every second
fes_parameters = {"model": DingModelFrequencyWithFatigue(stim_time=stim_time)}
ivp_parameters = {"final_time": 10}

ivp = IvpFes(fes_parameters, ivp_parameters)

result, time = ivp.integrate()

To personalize FES models to simulated or experimental force, Cocofest supports model identification using optimal control, through the OcpFesId class. Model parameters are treated as optimization parameters and identified by minimizing the difference between a tracked force (simulated with IvpFes, or your own experimental data) and the model's predicted force.

from cocofest import ModelMaker, OcpFesId, OcpFes
from cocofest.identification.identification_method import DataExtraction

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

# force_tracking: simulated (via IvpFes) or experimental (time, force) data to fit
force_at_node = DataExtraction.force_at_node_in_ocp(time, force, n_shooting, final_time)

x_bounds, x_init = OcpFesId.set_x_bounds(model=model, force_tracking=force_at_node)
u_bounds, u_init = OcpFesId.set_u_bounds(model=model, control_value=pulse_intensity_values, ...)

additional_key_settings = OcpFesId.set_default_values(model)
parameters, parameters_bounds, parameters_init = OcpFesId.set_parameters(
    parameter_to_identify=["a_rest", "km_rest", "tau1_rest", "tau2"],
    parameter_setting=additional_key_settings,
    use_sx=True,
)
OcpFesId.update_model_param(model, parameters)

ocp = OptimalControlProgram(bio_model=[model], x_bounds=x_bounds, x_init=x_init, u_bounds=u_bounds, u_init=u_init,
                             parameters=parameters, parameter_bounds=parameters_bounds, parameter_init=parameters_init, ...)
sol = ocp.solve()
identified_a_rest = sol.parameters["a_rest"][0]

Figure 3: Model identification optimization

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)

We are always looking for new contributors to help us improve Cocofest.
Feel free to check our contributing guidelines to get started, and please read our code of conduct beforehand.

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

The Cocofest companion paper is not published yet (submitted to JOSS).
Meanwhile, if you use Cocofest, please cite the software directly via its Zenodo archive: 10.5281/zenodo.17068808