******************************* Demo example - finger movements ******************************* The corresponding jupyter notebook script can be downloaded :download:`here <../scripts/demo_example_finger_movements.ipynb>`. The necessary input data (surfaces + muscles description) can be downloaded :download:`here: `. Import libraries ---------------- We start by importing necessary libraries. Myoelectric digital twin is part of the ``neurodec`` software. Apart from standard auxiliary libraries we also need ``meshio`` library to work with ``.stl`` surfaces .. code:: ipython3 %config Completer.use_jedi = False import glob import os import json import matplotlib.pyplot as plt import trimesh import numpy as np import neurodec as nd nd.mdt._API_URL = "http://35.173.188.106:8001/" nd.mdt.wait = False # Set to true to wait step by step. Load model surfaces ------------------- Subject’s forearm anatomy is represented by ``.stl`` meshes of the different tissue surfaces: upper skin, lower skin, individual muscles and bones. The surface filename for each muscle is stored in ``muscles_description.json`` file. We load each mesh with meshio package and create a ``nd.mdt.Surface`` for each tissue surface. .. code:: ipython3 def load_surface(filename: str): mesh = trimesh.load_mesh(filename) return mesh.vertices, mesh.faces MODEL_DIR = "/right-arm-template" muscles_description_fname = os.path.join(MODEL_DIR, 'muscles_description.json') with open(muscles_description_fname) as f: muscles_description = json.load(f) .. code:: ipython3 bones = [] for filename in glob.glob(os.path.join(MODEL_DIR, 'bones', '*.stl')): vertices, triangles = load_surface(filename) bones.append(nd.mdt.Surface.new(vertices, triangles, nd.mdt.SurfaceType.BONE)) .. code:: ipython3 muscles = [] for muscle in muscles_description['muscles']: fname = os.path.join(MODEL_DIR, 'muscles', muscle['fname']) vertices, triangles = load_surface(fname) muscles.append(nd.mdt.Surface.new(vertices, triangles, nd.mdt.SurfaceType.MUSCLE, label=muscle['label'])) .. code:: ipython3 skins = [] filename = os.path.join(MODEL_DIR, 'skin', 'inner.stl') vertices, triangles = load_surface(filename) skins.append(nd.mdt.Surface.new(vertices, triangles, nd.mdt.SurfaceType.INNER_SKIN)) filename = os.path.join(MODEL_DIR, 'skin', 'outer.stl') vertices, triangles = load_surface(filename) skins.append(nd.mdt.Surface.new(vertices, triangles, nd.mdt.SurfaceType.OUTER_SKIN)) Generate volume conductor ------------------------- Volume conductor is a function of anatomy surfaces and conductivities per tissue. Here, the default conductivity values are used. .. code:: ipython3 conductor = nd.mdt.Conductor.new(bones + muscles + skins) Generate bracelet electrodes ---------------------------- We generate a bracelet of electrodes made of ten rings. Each ring has 16 equidistant electrodes of 3mm radius. Distance between rings is 9 mm. Note that user can create electrodes by directly providing their centers and radii with ``nd.mdt.Electrode.new(location, radius)``. .. code:: ipython3 outer_skin_surface = skins[1] radius = 0.003 # in meters first_electrode_location = np.array([-0.247705, -0.083379, 0.838751]) rings_normal = np.array([0.17797643704463834, 0.19429176558661732, 0.9646632042750106]) n_rings = 10 distance_between_rings = 0.009 # in meters n_electrodes_per_ring = 16 bracelet = nd.mdt.ElectrodeBracelet.new( outer_skin_surface, radius, first_electrode_location, rings_normal, n_rings, distance_between_rings, n_electrodes_per_ring) You can save generated electrodes as spheres with ``trimesh`` for futher visualization .. code:: ipython3 print(bracelet.electrodes[0].location) for i, electrode in enumerate(bracelet.electrodes): sphere = trimesh.creation.icosphere() sphere.vertices = sphere.vertices * electrode.radius sphere.vertices = sphere.vertices + electrode.location sphere.export(MODEL_DIR + f'/elec_{i}.stl') Generate forward solver ----------------------- Forward solver is a function of volume conductor and electrodes .. code:: ipython3 forward_solution = nd.mdt.ForwardSolution.new(conductor, bracelet.electrodes) Generate fiber geometry ----------------------- User can directly provide the number of fibers, or it can be automatically estimated from the average fiber radius. Start and ending planes for fibers are stored in the ``muscle_description``. .. code:: ipython3 fibers_per_muscle = [] for i, muscle in enumerate(muscles): plane_origins = muscles_description['muscles'][i]['plane_origins'] plane_normals = muscles_description['muscles'][i]['plane_normals'] fibers = nd.mdt.Fibers.new(muscle, plane_origins, plane_normals) fibers_per_muscle.append(fibers) We also add individual finger “muscle” surfaces. Note, that those are not anatomical muscles but a subdivision of muscles controlling fingers so that individual finger motor units could be generated. .. code:: ipython3 finger_surfaces_description_fname = os.path.join(MODEL_DIR, 'individual_finger_surfaces.json') with open(finger_surfaces_description_fname) as f: finger_surfaces_description = json.load(f) for muscle in finger_surfaces_description['muscles']: # load extra surfaces for individual fingers fname = os.path.join(MODEL_DIR, 'muscles', muscle['fname']) vertices, triangles = load_surface(fname) muscles.append(nd.mdt.Surface.new(vertices, triangles, nd.mdt.SurfaceType.MUSCLE, label=muscle['label'])) # Generate corresponding fibers plane_origins = muscle['plane_origins'] plane_normals = muscle['plane_normals'] fibers = nd.mdt.Fibers.new(muscles[-1], plane_origins, plane_normals) fibers_per_muscle.append(fibers) Generate fiber basis -------------------- Fiber basis is a function of fibers and forward solution. .. code:: ipython3 fiber_basis_per_muscle = [] for i, fibers in enumerate(fibers_per_muscle): fiber_basis = nd.mdt.FiberBasis.new(fibers, forward_solution) fiber_basis_per_muscle.append(fiber_basis) Compute MUAPs ------------- To compute MUAPs for each muscle, fiber properties and motor units (MUs) should be generated and a sampling frequency (in Hz) should be provided. Here, fiber properties are generated using default tendon, velocity and neuromascular junction values. The number of MU should be specified. Here, default values for minimum and maximum relative MU areas are used. .. code:: ipython3 SAMPLING_FREQUENCY = 8000 n_motor_units = 150 muaps_per_muscle = [] for fibers, fiber_basis in zip(fibers_per_muscle, fiber_basis_per_muscle): fiber_properties = nd.mdt.FiberProperties.new(fibers) motor_units = nd.mdt.MotorUnits.new(fibers, n_motor_units) muaps = nd.mdt.MotorUnitsActionPotentials.new(fibers, fiber_basis, fiber_properties, motor_units, SAMPLING_FREQUENCY) muaps_per_muscle.append(muaps) Simulate finger movements ------------------------- Print the list of muscle labels .. code:: ipython3 for i, muscle in enumerate(muscles): print(i, muscle.label) .. parsed-literal:: 0 Extensor carpi radialis longus 1 Extensor carpi radialis brevis 2 Extensor carpi ulnaris 3 Flexor carpi radialis 4 Palmaris longus 5 Flexor carpi ulnaris humeral head 6 Flexor carpi ulnaris ulnar head 7 Pronator teres humeral head 8 Pronator teres ulnar head 9 Pronator quadratus 10 Supinator 11 Brachioradialis 12 Flexor digitorum superficialis humeroulnar head 13 Flexor digitorum superficialis radial head 14 Abductor pollicis longus 15 Anconeus 16 Extensor digiti minimi 17 Extensor digitorum 18 Extensor indicis 19 Extensor pollicis brevis 20 Extensor pollicis longus 21 Flexor digitorum profondus 22 Flexor pollicis longus 23 Extensor digitorum INDEX 24 Extensor digitorum MIDDLE 25 Extensor digitorum RING 26 Flexor digitorum profondus INDEX 27 Flexor digitorum profondus MIDDLE 28 Flexor digitorum profondus RING 29 Flexor digitorum profondus PINKY 30 Flexor digitorum superficialis INDEX 31 Flexor digitorum superficialis MIDDLE 32 Flexor digitorum superficialis RING 33 Flexor digitorum superficialis PINKY Select muscle indices for desired movements. Here are some examples: .. code:: ipython3 movements = { # Thumb flexion: FPL 'thumb_flex': [22], # Thumb abduction: APL, EPB 'thumb_abduct': [14, 19], # Thumb extension (lift): EPL 'thumb_ext': [20], # Index flexion: FDP index, FDS index 'index_flex': [26, 30], # Index extension: EI, EDC index 'index_ext': [18, 23], # Middle flexion: FDP middle, FDS middle 'middle_flex': [27, 31], # Middle extension: EDC middle 'middle_ext': [24], # Ring flexion: FDP ring, FDS ring 'ring_flex': [28, 32], # Ring extension: EDC ring 'ring_ext': [25], # Pinky flexion: FDP pinky, FDS pinky 'pinky_flex': [29, 33], # Pinky extension: EDM 'pinky_ext': [16] } Check the target muscles .. code:: ipython3 for movement, target_muscles in zip(movements.keys(), movements.values()): print(movement) for muscle_id in target_muscles: print('\t' + muscles[muscle_id].label) .. parsed-literal:: thumb_flex Flexor pollicis longus thumb_abduct Abductor pollicis longus Extensor pollicis brevis thumb_ext Extensor pollicis longus index_flex Flexor digitorum profondus INDEX Flexor digitorum superficialis INDEX index_ext Extensor indicis Extensor digitorum INDEX middle_flex Flexor digitorum profondus MIDDLE Flexor digitorum superficialis MIDDLE middle_ext Extensor digitorum MIDDLE ring_flex Flexor digitorum profondus RING Flexor digitorum superficialis RING ring_ext Extensor digitorum RING pinky_flex Flexor digitorum profondus PINKY Flexor digitorum superficialis PINKY pinky_ext Extensor digiti minimi Let’s define higher level function to simulate specific muscle activation patterns .. code:: ipython3 def assemble_emg(emg_list): max_length = np.max([e.shape[1] for e in emg_list]) emg_total = np.zeros([emg_list[0].shape[0], max_length]) for emg in emg_list: n_samples = emg.shape[1] emg_total[:, :n_samples] += emg return emg_total def single_muscle_activation(max_mvc, duration, peak_time, fs): '''Generates a single activation pattern smoothly rising to the max MVC and back to zero''' n_time_samples = int((peak_time + duration / 2.) * fs) t = np.linspace(-np.pi / 2, np.pi + np.pi / 2, int(duration * fs)) waveform = max_mvc * 0.5 * (np.sin(t) + 1.) i_0 = n_time_samples - len(waveform) muscle_activation = np.zeros(n_time_samples) muscle_activation[i_0:] = waveform return muscle_activation def generate_muscle_activations(n_muscles, mean_max_mvc, std_max_mvc, duration, peak_time, fs): '''Generates a multilple activation patterns with some variability in max MVC per muscle''' muscle_activations = [] for i in range(n_muscles): max_mvc = np.random.uniform(mean_max_mvc - std_max_mvc, mean_max_mvc + std_max_mvc) max_mvc = np.clip(max_mvc, 0., 100.) muscle_activations.append(single_muscle_activation(max_mvc, duration, peak_time, fs)) return muscle_activations def generate_movement(target_muscles, muaps_per_muscle, mean_max_mvc, std_max_mvc, duration, peak_time, fs, random_seed=42): '''Generates EMG corresponding to a specific movement''' # Generate muscle activations n_muscles = len(target_muscles) muscle_activations = generate_muscle_activations(n_muscles, mean_max_mvc, std_max_mvc, duration, peak_time, fs) # Generate impulse trains impulse_trains_per_muscle = [] for muscle_id, single_activation in zip(target_muscles, muscle_activations): impulse_trains_per_muscle.append(nd.mdt.ImpulseTrains.new(muaps_per_muscle[muscle_id], single_activation, random_seed=random_seed)) # Generate EMG for each muscle emg_per_muscle = [] for impulse_trains in impulse_trains_per_muscle: emg = nd.mdt.Electromyography.new(impulse_trains) emg.wait() # wait for the computation to finish before loading the EMG data emg_per_muscle.append(emg.data.copy()) # Sumup EMG signals of induvidual muscles emg_total = assemble_emg(emg_per_muscle) return muscle_activations, emg_total Example: 2 seconds long Index extension with 40%MVC peaking at 1.5 sec -------------------------------------------------------------------------- .. code:: ipython3 target_muscles = movements['index_ext'] mean_max_mvc = 40 # %MVC std_max_mvc = 10 # add some inter muscle variability in %MVC duration = 2. # sec peak_time = 1.5 # sec muscle_activations, emg = generate_movement(target_muscles, muaps_per_muscle, mean_max_mvc, std_max_mvc, duration, peak_time, SAMPLING_FREQUENCY, random_seed=42) .. code:: ipython3 # plot muscle activations plt.figure() for single_activation in muscle_activations: times = np.arange(len(single_activation)) / SAMPLING_FREQUENCY plt.plot(times, single_activation) # plot EMG times = np.arange(emg.shape[1]) / SAMPLING_FREQUENCY plot_shift = np.ptp(emg) rings_to_plot = [0, 4, 9] for ring in rings_to_plot: emg_to_plot = emg[ring*16:(ring+1)*16] plt.figure() plt.plot(times, emg_to_plot.T + plot_shift * np.arange(16)) plt.show() .. figure:: ../images/finger_movements/output_35_0.png Muscle activations .. figure:: ../images/finger_movements/output_35_1.png Electrodes ring 0 .. figure:: ../images/finger_movements/output_35_2.png Electrodes ring 4 .. figure:: ../images/finger_movements/output_35_3.png Electrodes ring 9 Example: Index extension followed by pinky extension ---------------------------------------------------- .. code:: ipython3 # Index target_muscles = movements['index_ext'] mean_max_mvc = 40 # %MVC std_max_mvc = 10 # add some inter muscle variability in %MVC duration = 2. # sec peak_time = 1.5 # sec _, emg_index = generate_movement(target_muscles, muaps_per_muscle, mean_max_mvc, std_max_mvc, duration, peak_time, SAMPLING_FREQUENCY, random_seed=42) # Pinky target_muscles = movements['pinky_ext'] mean_max_mvc = 40 # %MVC std_max_mvc = 10 # add some inter muscle variability in %MVC duration = 2. # sec peak_time = 3.5 # sec _, emg_pinky = generate_movement(target_muscles, muaps_per_muscle, mean_max_mvc, std_max_mvc, duration, peak_time, SAMPLING_FREQUENCY, random_seed=42) # Total emg_total = assemble_emg([emg_index, emg_pinky]) # plot EMG times = np.arange(emg_total.shape[1]) / SAMPLING_FREQUENCY plot_shift = np.ptp(emg_total) rings_to_plot = [0, 4, 9] for ring in rings_to_plot: emg_to_plot = emg_total[ring*16:(ring+1)*16] plt.figure() plt.plot(times, emg_to_plot.T + plot_shift * np.arange(16)) plt.show() .. figure:: ../images/finger_movements/output_37_0.png Electrodes ring 0 .. figure:: ../images/finger_movements/output_37_1.png Electrodes ring 4 .. figure:: ../images/finger_movements/output_37_2.png Electrodes ring 9 Example: Simultaneous index and thumb flexion --------------------------------------------- .. code:: ipython3 # Index target_muscles = movements['index_flex'] mean_max_mvc = 40 # %MVC std_max_mvc = 10 # add some inter muscle variability in %MVC duration = 2. # sec peak_time = 1.5 # sec _, emg_index = generate_movement(target_muscles, muaps_per_muscle, mean_max_mvc, std_max_mvc, duration, peak_time, SAMPLING_FREQUENCY, random_seed=42) # Thumb target_muscles = movements['thumb_flex'] mean_max_mvc = 40 # %MVC std_max_mvc = 10 # add some inter muscle variability in %MVC duration = 2. # sec peak_time = 1.5 # sec _, emg_thumb = generate_movement(target_muscles, muaps_per_muscle, mean_max_mvc, std_max_mvc, duration, peak_time, SAMPLING_FREQUENCY, random_seed=42) # Total emg_total = assemble_emg([emg_index, emg_thumb]) # plot EMG times = np.arange(emg_total.shape[1]) / SAMPLING_FREQUENCY plot_shift = np.ptp(emg_total) rings_to_plot = [0, 4, 9] for ring in rings_to_plot: emg_to_plot = emg_total[ring*16:(ring+1)*16] plt.figure() plt.plot(times, emg_to_plot.T + plot_shift * np.arange(16)) plt.show() .. figure:: ../images/finger_movements/output_39_0.png Electrodes ring 0 .. figure:: ../images/finger_movements/output_39_1.png Electrodes ring 4 .. figure:: ../images/finger_movements/output_39_2.png Electrodes ring 9