.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/06_simulate_force.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_06_simulate_force.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_06_simulate_force.py:


Force Generation
================

After generating **spike trains** from motor neuron pools, the next step is to simulate the **force output**
produced by the muscle. **MyoGen** provides a comprehensive force model based on the classic
Fuglevand et al. (1993) approach.

.. note::
    The **force model** converts spike trains into force by simulating individual motor unit twitches
    and their summation. Each motor unit has unique twitch characteristics (amplitude, duration)
    that depend on its recruitment threshold.

The force model includes:

* **Motor unit twitch characteristics**: Peak force and contraction time based on recruitment thresholds
* **Nonlinear gain modulation**: Force scaling based on inter-pulse intervals (discharge rate)
* **Temporal summation**: Realistic force buildup from overlapping twitches
* **Physiological scaling**: Realistic force ranges and dynamics

References
----------
.. [1] Fuglevand, A. J., Winter, D. A., & Patla, A. E. (1993).
    Models of recruitment and rate coding in motor-unit pools.
    Journal of Neurophysiology, 70(2), 782-797.

.. GENERATED FROM PYTHON SOURCE LINES 29-31

Import Libraries
----------------

.. GENERATED FROM PYTHON SOURCE LINES 31-46

.. code-block:: Python


    from pathlib import Path

    import numpy as np
    import matplotlib.pyplot as plt

    from myogen import simulator
    from myogen.simulator.core.force.force_model import ForceModel
    from myogen.utils.plotting.force import (
        plot_twitch_parameter_assignment,
        plot_twitches,
    )
    from myogen.utils.plotting import plot_spike_trains
    from myogen.utils.currents import create_trapezoid_current








.. GENERATED FROM PYTHON SOURCE LINES 47-59

Define Parameters
-----------------

The **force model** requires several key parameters:

- ``recruitment_thresholds``: Motor unit recruitment thresholds
- ``recording_frequency__Hz``: Sampling frequency for force simulation
- ``longest_duration_rise_time__ms``: Maximum twitch rise time
- ``contraction_time_range``: Range of contraction times across motor units

For demonstration, we'll also create input currents and simulate the complete
pipeline from current input to force output.

.. GENERATED FROM PYTHON SOURCE LINES 59-78

.. code-block:: Python


    # Create results directory
    save_path = Path("./results")
    save_path.mkdir(exist_ok=True)

    # Motor unit pool parameters
    n_motor_units = 100
    recruitment_range = 100  # Recruitment range (max_threshold / min_threshold)

    # Force model parameters
    recording_frequency__Hz = 10000.0  # 10 kHz sampling rate
    longest_duration_rise_time__ms = 90.0  # Maximum twitch rise time
    contraction_time_range = 3.0  # Contraction time range factor

    # Simulation parameters
    simulation_duration__ms = 10000.0  # 10 seconds
    timestep__ms = 0.1  # 0.1 ms time step
    t_points = int(simulation_duration__ms / timestep__ms)








.. GENERATED FROM PYTHON SOURCE LINES 79-85

Generate Recruitment Thresholds
--------------------------------

First, we need to generate recruitment thresholds for our motor unit pool.
These will determine both the activation order and the force characteristics
of each motor unit.

.. GENERATED FROM PYTHON SOURCE LINES 85-90

.. code-block:: Python


    recruitment_thresholds, _ = simulator.generate_mu_recruitment_thresholds(
        N=n_motor_units, recruitment_range__ratio=recruitment_range
    )








.. GENERATED FROM PYTHON SOURCE LINES 91-96

Create Force Model
------------------

The force model calculates individual motor unit twitch properties based on
recruitment thresholds and physiological scaling rules.

.. GENERATED FROM PYTHON SOURCE LINES 96-115

.. code-block:: Python


    force_model = ForceModel(
        recruitment_thresholds=recruitment_thresholds,
        recording_frequency__Hz=recording_frequency__Hz,
        longest_duration_rise_time__ms=longest_duration_rise_time__ms,
        contraction_time_range__unitless=contraction_time_range,
    )

    # Display force model statistics
    print(f"\nForce model statistics:")
    print(f"  - Number of motor units: {force_model._number_of_neurons}")
    print(f"  - Recruitment ratio: {force_model._recruitment_ratio:.1f}")
    print(
        f"  - Peak force range: {force_model.peak_twitch_forces__unitless[0]:.3f} - {force_model.peak_twitch_forces__unitless[-1]:.3f}"
    )
    print(
        f"  - Contraction time range: {force_model.contraction_times__samples[0]:.1f} - {force_model.contraction_times__samples[-1]:.1f} samples"
    )





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    Force model statistics:
      - Number of motor units: 100
      - Recruitment ratio: 100.0
      - Peak force range: 1.047 - 100.000
      - Contraction time range: 890.2 - 300.0 samples




.. GENERATED FROM PYTHON SOURCE LINES 116-122

Visualize Twitch Parameter Assignment
--------------------------------------

The force model assigns twitch parameters (peak force and contraction time)
to each motor unit based on its recruitment threshold. Let's visualize this
relationship for a subset of motor units.

.. GENERATED FROM PYTHON SOURCE LINES 122-146

.. code-block:: Python


    # Suppress font warnings to keep output clean
    import warnings
    import logging

    warnings.filterwarnings("ignore", message=".*Font family.*not found.*")
    warnings.filterwarnings("ignore", message=".*findfont.*")
    logging.getLogger("matplotlib.font_manager").setLevel(logging.ERROR)

    plt.figure(figsize=(8, 12))

    ax1 = plt.subplot(2, 1, 1)
    plot_twitch_parameter_assignment(
        force_model, ax1, [1, 50, 100], flip_x=True, apply_default_formatting=True
    )
    ax1.set_title("Twitch Parameter Assignment")

    ax2 = plt.subplot(2, 1, 2)
    plot_twitches(force_model, ax2, apply_default_formatting=True)
    ax2.set_title("Motor Unit Twitches")

    plt.tight_layout()
    plt.show()




.. image-sg:: /auto_examples/images/sphx_glr_06_simulate_force_001.png
   :alt: Twitch Parameter Assignment, Motor Unit Twitches
   :srcset: /auto_examples/images/sphx_glr_06_simulate_force_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 147-152

Generate Input Current
----------------------

To demonstrate force generation, we'll create a trapezoid current that will
drive motor unit recruitment and firing patterns.

.. GENERATED FROM PYTHON SOURCE LINES 152-187

.. code-block:: Python


    # Parameters for trapezoid current
    trap_amplitude = 100.0  # Peak amplitude
    trap_rise_time = 1000.0  # Rise duration (ms)
    trap_plateau_time = 3000.0  # Plateau duration (ms)
    trap_fall_time = 1000.0  # Fall duration (ms)
    trap_offset = 0.0  # Baseline current
    trap_delay = 2500.0  # Initial delay (ms)

    trapezoid_current = create_trapezoid_current(
        n_pools=1,
        t_points=t_points,
        timestep__ms=timestep__ms,
        amplitudes__muV=[trap_amplitude],
        rise_times__ms=[trap_rise_time],
        plateau_times__ms=[trap_plateau_time],
        fall_times__ms=[trap_fall_time],
        offsets__muV=[trap_offset],
        delays__ms=[trap_delay],
    )

    # Plot the input current
    plt.figure(figsize=(12, 4))
    with plt.xkcd():
        time_vector = np.arange(t_points) * timestep__ms / 1000  # Convert to seconds
        plt.plot(time_vector, trapezoid_current[0], "g-", linewidth=2)
        plt.xlabel("Time (s)")
        plt.ylabel("Current (µV)")
        plt.title("Input Current Waveform")
        plt.grid(True, alpha=0.3)
        plt.xlim(0, simulation_duration__ms / 1000)

    plt.tight_layout()
    plt.show()




.. image-sg:: /auto_examples/images/sphx_glr_06_simulate_force_002.png
   :alt: Input Current Waveform
   :srcset: /auto_examples/images/sphx_glr_06_simulate_force_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 188-193

Create Motor Neuron Pool and Generate Spike Trains
---------------------------------------------------

Now we'll create a motor neuron pool and generate spike trains in response
to the input current.

.. GENERATED FROM PYTHON SOURCE LINES 193-207

.. code-block:: Python


    # Create motor neuron pool
    motor_neuron_pool = simulator.MotorNeuronPool(recruitment_thresholds)

    # Generate spike trains
    spike_trains_matrix, active_neuron_indices, data = (
        motor_neuron_pool.generate_spike_trains(
            input_current__matrix=trapezoid_current,
            timestep__ms=timestep__ms,
            noise_mean__nA=0.0,
            noise_stdev__nA=0.0,
        )
    )





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    Creating motor neuron pools:   0%|          | 0/1 [00:00<?, ?pool/s]
    Creating motor neuron pools: 100%|██████████| 1/1 [00:00<00:00, 24.62pool/s]

    Injecting currents into pools:   0%|          | 0/1 [00:00<?, ?pool/s]

    Adding noise to neurons in pool:   0%|          | 0/100 [00:00<?, ?neuron/s]

    Adding noise to neurons in pool:  21%|██        | 21/100 [00:00<00:00, 204.81neuron/s]

    Adding noise to neurons in pool:  42%|████▏     | 42/100 [00:00<00:00, 203.57neuron/s]

    Adding noise to neurons in pool:  63%|██████▎   | 63/100 [00:00<00:00, 204.59neuron/s]

    Adding noise to neurons in pool:  84%|████████▍ | 84/100 [00:00<00:00, 205.31neuron/s]

                                                                                          
    Injecting currents into pools: 100%|██████████| 1/1 [00:00<00:00,  1.63pool/s]
    Injecting currents into pools: 100%|██████████| 1/1 [00:00<00:00,  1.63pool/s]




.. GENERATED FROM PYTHON SOURCE LINES 208-213

Visualize Spike Trains
-----------------------

Let's visualize the generated spike trains to see the motor unit recruitment
and firing patterns.

.. GENERATED FROM PYTHON SOURCE LINES 213-226

.. code-block:: Python


    plt.figure(figsize=(12, 8))

    with plt.xkcd():
        ax = plt.gca()
        plot_spike_trains(spike_trains_matrix, timestep__ms, [ax])
        ax.set_title("Motor Unit Spike Trains")
        ax.set_xlabel("Time (s)")
        ax.set_ylabel("Motor Unit Number")

    plt.tight_layout()
    plt.show()




.. image-sg:: /auto_examples/images/sphx_glr_06_simulate_force_003.png
   :alt: Motor Unit Spike Trains
   :srcset: /auto_examples/images/sphx_glr_06_simulate_force_003.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 227-233

Generate Force Output
---------------------

Now we can use the force model to convert the spike trains into force output.
The model will simulate individual motor unit twitches and their temporal
summation to produce the total muscle force.

.. GENERATED FROM PYTHON SOURCE LINES 233-243

.. code-block:: Python


    # Generate force from spike trains
    force_output = force_model.generate_force(spike_trains_matrix)

    # Add some realistic noise to the force signal
    noise_level = 0.015  # 0.5% of mean force
    noisy_force = force_output[0] + np.random.randn(
        len(force_output[0])
    ) * noise_level * np.mean(force_output[0])





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


     Twitch trains are generated:   0%|          | 0/100 [00:00<?, ?it/s]
     Twitch trains are generated:   9%|▉         | 9/100 [00:00<00:01, 85.19it/s]
     Twitch trains are generated:  18%|█▊        | 18/100 [00:00<00:00, 85.47it/s]
     Twitch trains are generated:  27%|██▋       | 27/100 [00:00<00:00, 85.71it/s]
     Twitch trains are generated:  36%|███▌      | 36/100 [00:00<00:00, 86.43it/s]
     Twitch trains are generated:  45%|████▌     | 45/100 [00:00<00:00, 87.06it/s]
     Twitch trains are generated:  55%|█████▌    | 55/100 [00:00<00:00, 88.14it/s]
     Twitch trains are generated:  65%|██████▌   | 65/100 [00:00<00:00, 89.25it/s]
     Twitch trains are generated:  75%|███████▌  | 75/100 [00:00<00:00, 91.45it/s]
     Twitch trains are generated:  85%|████████▌ | 85/100 [00:00<00:00, 92.30it/s]
     Twitch trains are generated:  95%|█████████▌| 95/100 [00:01<00:00, 92.45it/s]
     Twitch trains are generated: 100%|██████████| 100/100 [00:01<00:00, 89.92it/s]




.. GENERATED FROM PYTHON SOURCE LINES 244-249

Visualize Force Output
----------------------

Let's plot the generated force alongside the input current to see how
the muscle responds to the electrical stimulation.

.. GENERATED FROM PYTHON SOURCE LINES 249-287

.. code-block:: Python


    plt.figure(figsize=(12, 10))

    with plt.xkcd():
        # Plot input current
        ax1 = plt.subplot(3, 1, 1)
        ax1.plot(
            time_vector, trapezoid_current[0], "g-", linewidth=2, label="Input Current"
        )
        ax1.set_ylabel("Current (µV)")
        ax1.set_title("Input Current and Force Response")
        ax1.grid(True, alpha=0.3)
        ax1.legend()
        ax1.set_xlim(0, simulation_duration__ms / 1000)

        # Plot clean force
        ax2 = plt.subplot(3, 1, 2)
        ax2.plot(time_vector, force_output[0], "b-", linewidth=2, label="Clean Force")
        ax2.set_ylabel("Force (N)")
        ax2.set_title("Simulated Force Output")
        ax2.grid(True, alpha=0.3)
        ax2.legend()
        ax2.set_xlim(0, simulation_duration__ms / 1000)

        # Plot noisy force (more realistic)
        ax3 = plt.subplot(3, 1, 3)
        ax3.plot(
            time_vector, noisy_force, "r-", linewidth=1, alpha=0.8, label="Noisy Force"
        )
        ax3.set_xlabel("Time (s)")
        ax3.set_ylabel("Force (N)")
        ax3.set_title("Realistic Force Output (with noise)")
        ax3.grid(True, alpha=0.3)
        ax3.legend()
        ax3.set_xlim(0, simulation_duration__ms / 1000)

    plt.tight_layout()
    plt.show()



.. image-sg:: /auto_examples/images/sphx_glr_06_simulate_force_004.png
   :alt: Input Current and Force Response, Simulated Force Output, Realistic Force Output (with noise)
   :srcset: /auto_examples/images/sphx_glr_06_simulate_force_004.png
   :class: sphx-glr-single-img






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (1 minutes 1.404 seconds)


.. _sphx_glr_download_auto_examples_06_simulate_force.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: 06_simulate_force.ipynb <06_simulate_force.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: 06_simulate_force.py <06_simulate_force.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: 06_simulate_force.zip <06_simulate_force.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_