Open In Colab

Kinematics and Kinetics of Walking

We can’t understand anything about how the nervous system organizes and controls a movement until we can understand the movement itself and the biomechanics of the body controlled by the nervous system.

Kinematics is the study and quantification of how objects move. In working through the contents of this notebook, you will learn more about movement and how it is measured and analyzed. You will also practice using data processing and visualization techniques common to neuroscience, physics, biology, and many data science professions.

Specifically, you will learn about how (some) human bodies walk in the context of current standard practices in the fields of biomechanics. The dataset you will explore was published by Fukuchi et al (2018). They placed markers placed at specific locations along the body were tracked during walking (either overland or on a treadmill). Note that the nomenclature for the marker points in this dataset vary slightly from the marker placement resource referenced in Fukuchi et al (2018). See marker key in slides from class. Force plates were used to quantify ground contact forces during walking (“kinetics”).

https://github.com/neurologic/Neurobiology-Movement/blob/main/images/kinematics_markers_Fukuchi.png?raw=true

Your work on this notebook will scaffold discussions throughout the course about kinematics of a variety of movements, the challenges faced by the nervous system in controlling movement, and emerging technologies that are revolutionizing ways of quantifying movement.

Setup

Complete the following preparatory steps to set up the Colab notebook environment for your work.

#@title {display-mode: "form"}

#@markdown RUN this cell to set up the notebook (import packages, define functions, etc)
# In Python, anything with a "#" in front of it is code annotation,
# and is not read by the computer.
# You can run a cell (this box) by pressing ctrl-enter or shift-enter.
# You can also run a cell by clicking the play button in the menu bar 
# at the top of the page (single right arrow, not double).
# Click in this cell and then press shift and enter simultaneously.
# This print function below allows us to generate a message.
print('Nice work!')

# No need to edit anything in this code cell
#################################

import numpy as np
import pandas as pd
from scipy import ndimage
from scipy.signal import find_peaks
from copy import deepcopy
import math
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans
from pathlib import Path
import matplotlib.pyplot as plt
import plotly.express as px
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from ipywidgets import widgets,interactive, interactive_output
import matplotlib.gridspec as gridspec
from IPython.display import display
import csv
from scipy.signal import hilbert,medfilt,resample
from scipy.io import wavfile
from sklearn.decomposition import PCA
import scipy
import seaborn as sns
from datetime import datetime,timezone,timedelta
pal = sns.color_palette(n_colors=15)
pal = pal.as_hex()

%config InlineBackend.figure_format = 'retina'
plt.style.use("https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/nma.mplstyle")

# datafolder = "/content/drive/Shareddrives/BIOL358/Data/WBDSascii"
# datafolder = "/content/drive/MyDrive/Classroom/BIOL358: Motor Systems/Data/WBDSascii"
# datafolder = '/Users/kperks/OneDrive - wesleyan.edu/Teaching/MotorSystems_SP22/Data/WBDSascii/'

fs_grf = 300
fs_mkr = 150

def cart2pol(x, y):
    rho = np.sqrt(x**2 + y**2)
    phi = np.arctan2(y, x)
    return(rho, phi)
    
print('Task completed at ' + str(datetime.now(timezone(-timedelta(hours=5)))))
#@title {display-mode: "form"}

#@markdown RUN this cell to download and unzip the raw data from "figshare".
#@markdown Note that this step may take a few minutes. If the "run cell" button is still spinning, then it is still working.

!wget https://figshare.com/ndownloader/files/10058986

# use --no-check-certificate to ignore SSL certificate errors that might occur on websites with an expired SSL certificate

!unzip -q /content/10058986 -d /content/WBDSascii 
datafolder = "/content/WBDSascii"

# !unzip -q ./10058986 -d ./WBDSascii 
# datafolder = "./WBDSascii"

print('Task completed at ' + str(datetime.now(timezone(-timedelta(hours=5)))))

The dataset should now be located in your Colab file directory in a folder called “WBDSascii”. Check that it is there before proceeding.

Dataset Selection

Run the code cell below to activate dropdown menus for selecting which person, environment, and speed you want to analyze (note different speed selection options based on the environment chosen).

NOTE: You only need to activate these menus once per active colab session.

#@title {display-mode: "form"}

#@markdown Activate Dropdown Menus if needed.
 
w_subj = widgets.Select(
    options=['select person (01-42)'] + [format(x,'02d') for x in range(1,43,1)],
    value='select person (01-42)',
    # rows=10,
    disabled=False
)
w_env = widgets.Select(
    options=['select environment','Overland','Treadmill'],
    value='select environment',
    rows=3,disabled=False
)
w_Ospeed = widgets.Select(
    options=['select speed for overland','Comfortable','Fast','Slow'],
    value='select speed for overland',
    rows=4,disabled=False
)
w_Tspeed = widgets.Select(
    options=['select speed for treadmill (01-08)'] + [format(x,'02d') for x in range(1,9,1)],
    value='select speed for treadmill (01-08)',
    rows=5,disabled=False
)
display(w_subj,w_env,w_Ospeed,w_Tspeed)
print('Dropdowns generated at ' + str(datetime.now(timezone(-timedelta(hours=5)))))

Load Data

#@title {display-mode: "form"} 

#@markdown Run this code cell to load the data you have selected.

if w_env.value=='Overland':
    speed = w_Ospeed.value[0]
    trial = '01'
    filename = 'WBDS' + str(w_subj.value) + 'walk' + w_env.value[0] + trial + speed 

if w_env.value=='Treadmill':
    speed = w_Tspeed.value
    filename = 'WBDS' + str(w_subj.value) + 'walk' + w_env.value[0] + speed 


filepath = Path(datafolder) / (filename + 'mkr.txt')
if filepath.exists():
    dataframe = pd.read_csv(filepath, sep='\t')
    print('Kinematic Data loaded at ' + str(datetime.now(timezone(-timedelta(hours=5)))))

if filepath.exists()==False:
    print('Data selection does not exist. Check selections and correct errors.')
    print('Task not yet completed, but cell ran at ' + str(datetime.now(timezone(-timedelta(hours=5)))))    

filepath_grf = Path(datafolder) / (filename + 'grf.txt')
if filepath_grf.exists():
    grf_trial = pd.read_csv(filepath_grf,sep='\t')
    grf_trial['Time']/=fs_grf
    print('Kinetic Data loaded at ' + str(datetime.now(timezone(-timedelta(hours=5)))))

if filepath_grf.exists()==False:
    print('Data selection does not exist. Check selections and correct errors.')
    print('Task not yet completed, but cell ran at ' + str(datetime.now(timezone(-timedelta(hours=5)))))
    
filename_norm = 'WBDS' + str(w_subj.value) + 'walk' + w_env.value[0] + speed 

filepath_norm = Path(datafolder) / (filename_norm + 'ang.txt')
ang_norm = pd.read_csv(filepath_norm,sep = '\t')
print('Processed kinematics data loaded at ' + str(datetime.now(timezone(-timedelta(hours=5)))))


filepath_norm = Path(datafolder) / (filename_norm + 'knt.txt')
grf_norm = pd.read_csv(filepath_norm,sep = '\t')
print('Processed ground fource data loaded at ' + str(datetime.now(timezone(-timedelta(hours=5)))))

RAW data

Kinematics

You can select different kinematics variables and re-make the following plot as many times as you want without re-loading the dataset.

Follow your curiosities about what kinematics are related to each other and in what ways. There are many to choose from. Explore both ‘ipsilateral’ and ‘contralateral’ comparisons.

#@title {display-mode: "form"}

#@markdown Run this code cell to get a list of kinematic variables to choose from.

w_markers = widgets.SelectMultiple(
    options=list(dataframe.columns)[1:],
    rows=15)

print('')
print('')
print('Select which kinematics variables you want to plot from the following list:')
display(w_markers)
print('')

print('Dropdown menus created at ' + str(datetime.now(timezone(-timedelta(hours=5)))))
#@markdown Run this code cell to plot the selected kinematics variables


fig = go.Figure()
if len(w_markers.value)>0:
    for marker in list(w_markers.value):
        fig.add_trace(go.Scatter(x = dataframe['Time'], y = dataframe[marker],name=marker))
    print('Task completed at ' + str(datetime.now(timezone(-timedelta(hours=5)))))

if len(w_markers.value)==0:
    print('Need to select markers to plot')
    print('Task not completed yet, but cell run at ' + str(datetime.now(timezone(-timedelta(hours=5)))))
    
fig.update_layout(xaxis_title="time(seconds)", yaxis_title='position (centimeters)',width=800, height=500)

Kinetics

Select a subset of ground force measurements to examine.

Treadmill environments have more ground force data than overland.

#@title {display-mode: "form"}

#@markdown Run this code cell to get a list of kinetic variables to choose from.

w_fp = widgets.SelectMultiple(
    options=list(grf_trial.columns)[1:],
    rows=10)

print('')
print('Select which kinetic variables you want to plot from the following list:')
display(w_fp)
print('Dropdown menus created at ' + str(datetime.now(timezone(-timedelta(hours=5)))))
#@title {display-mode: "form"}

#@markdown Run this code cell to plot the kinetics data you selected.

fig = go.Figure()
if len(w_fp.value)>0:
    for fp in list(w_fp.value):
        fig.add_trace(go.Scatter(x = grf_trial['Time'], y = grf_trial[fp],name=fp))
    print('Plot created at ' + str(datetime.now(timezone(-timedelta(hours=5)))))

if len(w_fp.value)==0:
    print('Need to select markers to plot')
    print('Task not completed, but cell run at ' + str(datetime.now(timezone(-timedelta(hours=5)))))

fig.update_layout(xaxis_title="time(seconds)", yaxis_title='force',width=800, height=500)

Compare across trials

First, load data from a treadmill condition.

The treadmill condition has more trials.

Each “trial” is be defined as a single gait cycle. The onset of a gait cycle is determined based on kinetic data (ground reaction forces measured by the force plates).

Based on the ground force measurements you just plotted in the last code cell, use your mouse hover to determine the start times of 3 trials. In the code cell below, use the form to create a list of trial start times and an approximate gait cycle period (seconds).

The formatting of the list must be trial start times (in seconds) separated by commas. The whole list must be surrounded by hard brackets. For example: [1.24,5.67,8.98]

After filling in the necessary information, run the code cell to overlay the kinetic data across those trials in a single plot.

#@title {display-mode: "form"}

#@markdown Type a list of trial times below. 
trials = [None] #@param 

#@markdown Type the gait cycle period below (in seconds). 
period =  None #@param 

#@markdown After you have created a list of trial start times and specified a period,
#@markdown run this cell to make a plot of the overlaid kinetics data. 

offset = 0.5
persamps = int(period*fs_grf)+int(offset*fs_grf)
trial_inds = [int((t-offset)*fs_grf) for t in trials]
xtime = np.linspace(-offset,period,persamps)


fig = go.Figure()

if len(w_fp.value)>0:
    for fp in list(w_fp.value):
      for t in trial_inds:
        fig.add_trace(go.Scatter(x = xtime, y = grf_trial[fp][t:t+persamps],name=fp+' at trial ' +str(round(t/fs_grf,2)+offset)))
    print('Task completed at ' + str(datetime.now(timezone(-timedelta(hours=5)))))

if len(w_fp.value)==0:
    print('Need to select markers to plot')
    print('Task not completed, but cell run at ' + str(datetime.now(timezone(-timedelta(hours=5)))))

fig.update_layout(xaxis_title="time(seconds)", yaxis_title='force',width=800, height=500)



# #@markdown Run this code cell to plot the kinematics measurements on each trial overlaid.

# offset = 0.5
# persamps = int(period*fs_mkr)+int(offset*fs_mkr)
# trial_inds = [int((t-offset)*fs_mkr) for t in trials]
# xtime = np.linspace(-offset,period,persamps)

# fig = go.Figure()

# if len(w_markers.value)>0:
#     for marker in list(w_markers.value):
#       for t in trial_inds:
#         fig.add_trace(go.Scatter(x = xtime, y = dataframe[marker][t:t+persamps],name=str(marker)+' at trial ' +str(round(t/fs_mkr,2))))
#     print('Task completed at ' + str(datetime.now(timezone(-timedelta(hours=5)))))

# if len(w_markers.value)==0:
#     print('Need to select markers to plot')
#     print('Task not completed, but cell run at ' + str(datetime.now(timezone(-timedelta(hours=5)))))

# fig.update_layout(xaxis_title="time(seconds)", yaxis_title='position',width=800, height=500)

What do you notice across trials?

Time-normalized (processed) data

Now that you have a sense of the data, let’s look at trial-averaged and time-normalized data. The data is “time-normalized” so that the kinematics and kinetics are a function of the gait cycle (from start to stop - 0 to 100% of the gait cycle) instead of time.

Think about why time normalization important for an analysis of kinematics and kinetics.

Kinematics and Kinetics

#@title {display-mode: "form"} 

#@markdown Run this cell to explore the average kinematics and kinetics per gait cycle. 

select_kinematics = widgets.SelectMultiple(
    options=list(ang_norm.columns)[1:], # start with a single trial on a single bout... it will update when runs ; old: np.arange(len(trial_times)),
    value=[list(ang_norm.columns)[1:][0]],
    rows=10,
    description='Kinematics',
    style = {'description_width': '100px'},
    disabled=False
)

select_kinetics = widgets.SelectMultiple(
    options=list(grf_norm.columns)[1:], # start with a single trial on a single bout... it will update when runs ; old: np.arange(len(trial_times)),
    value=[list(grf_norm.columns)[1:][0]],
    rows=10,
    description='Kinetics',
    style = {'description_width': '100px'},
    disabled=False
)

ui_selections = widgets.HBox([select_kinematics, select_kinetics])

def update_plot(kinematics_,kinetics_):
    fig, ax = plt.subplots(figsize=(10,7),num=1,nrows=2,ncols=1); #specify figure number so that it does not keep creating new ones
    
    for m_ in kinematics_:
        ax[0].plot(ang_norm['Time'],ang_norm[m_],label=m_)

    ax[0].legend(bbox_to_anchor=(1.05, 1))
    ax[0].set_ylabel('degrees')
        
    for n_ in kinetics_:
        ax[1].plot(grf_norm['Time'],grf_norm[n_],label=n_)

    ax[1].legend(bbox_to_anchor=(1.05, 1))
    ax[1].set_ylabel('newton-meters')
    ax[1].set_xlabel('percent gait cycle')

# w_normalized_plots_ = interactive(update_plot, kinematics_=select_kinematics, 
#                         kinetics_=select_kinetics);
# display(w_normalized_plots_)

w = interactive_output(update_plot, {'kinematics_':select_kinematics,
                                     'kinetics_':select_kinetics
                                     });
display(ui_selections,w)

Ground Reaction Force

#@title {display-mode: "form"} 

#@markdown Run this cell to plot the ground reaction force vector direction and magnitude


fig = plt.figure(figsize=(15,10),num=1)

gs = gridspec.GridSpec(2, 2, figure=fig)
ax1 = fig.add_subplot(gs[0, 0])
ax2 = fig.add_subplot(gs[1, 0])
ax3 = fig.add_subplot(gs[:, 1], polar=True) # identical to ax3 = plt.subplot(gs.new_subplotspec((0, 0), rowspan=2))


rho,phi = cart2pol(grf_norm['RGRFX'],grf_norm['RGRFY'])
ax1.scatter(grf_norm['Time'],np.degrees(phi),label='right',s=10,color='blue')
ax2.scatter(grf_norm['Time'],rho,label='right',s=10,color='blue')
ax3.scatter(phi,rho, color='blue', alpha=0.5,s=5)

rho,phi = cart2pol(grf_norm['LGRFX'],grf_norm['LGRFY'])
ax1.scatter(grf_norm['Time'],np.degrees(phi),label='left',s=10,color='red')
ax2.scatter(grf_norm['Time'],rho,label='left',s=10,color='red')
ax3.scatter(phi,rho, color='red', alpha=0.5,s=5)

ax1.legend();
ax1.set_ylabel('angle')
ax2.set_ylabel('magnitude')
ax2.set_xlabel('percent gait cycle');


NICE WORK!



For users not in this course: This notebook assumes that you have access to the BIOL358 shared class drive in which the dataset from Fukuchi et al (2018) has been stored. To change the filepath, or aquire the dataset via direct download (if that is possible), edit the first code cell in Step #2 (“TASK: RUN this cell to set up the notebook (import packages, etc)” and change the datafolder path variable. And/or edit the code cells in which the data is imported.