Experimental SSCs

The analysis of this page corresponds to the sections ‘stretch-shortening cycles’.

Because the analysis is relatively extensive, it is split into three separate scripts included on this page:

1. Compute AMPO for each cycle for every rat, condition, and trial. The results are stored in a spreadsheet.
2. Perform statistical analyses on the effects of cycle frequency and FTS on the measured AMPO.
3. Calculate differences in measured AMPO between conditions.

Custom functions used:

• stimulation.get_stim_timing(time, stim)

  Detects stimulation pulse trains in a signal and returns their onset and offset time.

• stats.pdiff(x,y)

  Compute the percentage difference with x w.r.t. y.

Make table

"""
This script computes experimentally measured AMPO for every stimulated SSC
cycle and stores the values in the per-specimen spreadsheet.

For each rat, condition and trial, stimulation timing is used to identify the
cycles. AMPO is computed from the work-loop area of each stimulated cycle.
Saving is disabled by default and controlled with `do_save`.
"""

#%% Load packages & set directories
import os, sys, scipy, openpyxl
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path

# Set directories
cwd = Path.cwd()
baseDir = cwd.parent.parent
dataDir = baseDir / 'data'
funcDir = baseDir / 'analysis' / 'functions'
sys.path.append(str(funcDir))

import stimulation

plt.close('all')

#%% Compute AMPO per cycle and store in sheet
do_save = False # set to true to save
for mus in ['GMe1', 'GMe2', 'GMe3']:
    AMPOarray = np.nan*np.empty((12,13))
    for exp in ['SSC_PA', 'SSC_PB']:
        for iCond in range(1,14):
            if exp == 'SSC_PA':
                cf = [1,2,3,4,5,3,3,3,3,5,4,2,1][iCond-1]
            elif exp == 'SSC_PB':
                cf = [1,1.5,2,2.5,3,2,2,2,2,3,2.5,1.5,1][iCond-1]
            for iTrial in [1,2]:
                if exp == 'SSC_PA':
                    iStartRow = [0,3][iTrial-1]
                elif exp == 'SSC_PB':
                    iStartRow = [6,9][iTrial-1]
                    
                filepath = os.path.join(dataDir,mus,'dataExp',exp,f'{mus}_{exp}{iCond:02d}_{iTrial}.csv')
                try:
                    data = pd.read_csv(filepath).T.to_numpy()
                    time, lmtc, stim, fsee, *_ = data
                    t_stimOn, t_stimOff = stimulation.get_stim_timing(time, stim)
                    
                    # Determine samples per cycle do this as follows:
                    # 1) Determine time between stimulation onsets and stimulation offsets
                    tCycle = np.mean(np.diff(t_stimOn))/2 + np.mean(np.diff(t_stimOff))/2
                    # 2) From time to number of samples
                    nCycle = int(tCycle/np.mean(np.diff(time)))
                    
                    # Determine iMax:
                    iMax = scipy.signal.find_peaks(lmtc,distance=nCycle*0.95)[0]
                    
                    if len(iMax) != 6:
                        if mus == 'GMe3' and exp=='SSC_PB' and iCond==1 and iTrial==2:
                            pass # checked: iMax @ end missing 
                        elif mus == 'GMe3' and exp=='SSC_PB' and iCond==3 and iTrial==1:
                            pass # checked: iMax @ end missing 
                        elif mus == 'GMe3' and exp=='SSC_PB' and iCond==5 and iTrial==1:
                            pass # checked: iMax @ end missing 
                        elif mus == 'GMe3' and exp=='SSC_PB' and iCond==7 and iTrial==1:
                            pass # checked: iMax @ end missing 
                        else:
                            breakpoint()
                    
                    nSamples = np.diff(iMax)[1:4]
                    if abs(nSamples-2000/cf).max() > 2:
                        breakpoint() # something goes wrong with findpeaks!
                    
                    for i in range(1,4):
                        iSel = slice(iMax[i],iMax[i+1])
                        AMPO = -scipy.integrate.trapezoid(fsee[iSel], lmtc[iSel])/(time[iSel][-1]-time[iSel][0])
                        if AMPO < 0:
                            # This should only print:
                                # GMe1, SSC_PB, iCond=2, iTrial=2, iCycle=3
                                # For some reason the 3rd cycle did not get stimulation here..
                            print(f'AMPO<0 for mus={mus}, exp={exp}, iCond={iCond}, iTrial={iTrial}, iCycle={i}')
                        AMPOarray[iStartRow+i-1,iCond-1] = AMPO*1e3 # to mW
                except Exception:
                    pass
    
    df = pd.DataFrame(AMPOarray)
    
    if do_save is True:
        # Path to existing Excel file
        filepath = os.path.join(dataDir,mus,f'{mus}_dataAMPO.xlsx')
        
        # Open the existing Excel file using openpyxl
        wb = openpyxl.load_workbook(filepath)
        
        # Select the sheet (you can also use wb[sheet_name] if you know the name)
        ws = wb.active  # or wb['Sheet1'] to select a specific sheet by name
        
        # Step 1: Define the starting position (row, column)
        start_row = 2  # Starting row where you want to insert data
        start_col = 6  # Starting column where you want to insert data (1 = 'A', 2 = 'B', etc.)
        
        # Step 2: Write the DataFrame to specific rows and columns
        for r_idx, row in enumerate(df.itertuples(index=False), start=start_row):
            for c_idx, value in enumerate(row, start=start_col):
                ws.cell(row=r_idx, column=c_idx, value=value)
        
        # Step 3: Save the modified workbook
        wb.save(filepath)
        wb.close()

Statistics

"""
This script tests whether cycle frequency and FTS affect experimentally
measured AMPO.

The AMPO spreadsheet is reshaped into a long-format dataframe for the four
subsets shown in the experimental AMPO figure. For each subset, a linear
mixed-effects model with linear and quadratic terms is fitted, using specimen
as random effect.
"""

#%% Load packages & set directories
import os, sys
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.formula.api import mixedlm
from pathlib import Path

# Set directories
cwd = Path.cwd()
baseDir = cwd.parent.parent
dataDir = baseDir / 'data'
funcDir = baseDir / 'analysis' / 'functions'
sys.path.append(str(funcDir))

#%% Make dataframe
data = {
    'specimen': [],
    'cf': [],
    'fts': [],
    'mle': [],
    'trial': [],
    'AMPO': [],
    'subset': [],
}
    
for iMus,mus in enumerate(['GMe1','GMe2','GMe3']):
    dataDirMus = os.path.join(dataDir,mus,'')
    dataExp = pd.read_excel(dataDirMus+str(mus)+'_dataAMPO.xlsx').to_numpy()[:,5:].astype(float)
    dataExp[dataExp < 1] = np.nan
        
    # Panel A: 4mm MLE - effect of CF
    cf = np.matlib.repmat([1,2,3,4,5],6,1) # cf
    fts = np.matlib.repmat([0.5],6,5) # fts
    mle = np.matlib.repmat([4],6,5)
    trial = np.repeat([1, 2], 3)[:, None].repeat(5, axis=1)
    AMPO = dataExp[0:6,0:5]
    
    data['specimen'].extend([mus]*30)
    data['cf'].extend(cf.flatten())
    data['fts'].extend(fts.flatten())
    data['mle'].extend(mle.flatten())
    data['trial'].extend(trial.flatten())
    data['AMPO'].extend(AMPO.flatten())
    data['subset'].extend(['A']*30)
    
    # Panel B: 8mm MLE - effect of FTS
    cf = np.matlib.repmat([1,1.5,2,2.5,3],6,1) # cf
    fts = np.matlib.repmat([0.5],6,5) # fts
    mle = np.matlib.repmat([8],6,5)
    trial = np.repeat([1, 2], 3)[:, None].repeat(5, axis=1)
    AMPO = dataExp[6:12,0:5]
    
    data['specimen'].extend([mus]*30)
    data['cf'].extend(cf.flatten())
    data['fts'].extend(fts.flatten())
    data['mle'].extend(mle.flatten())
    data['trial'].extend(trial.flatten())
    data['AMPO'].extend(AMPO.flatten())
    data['subset'].extend(['B']*30)
    
    # Panel C: 4mm MLE - effect of CF
    cf = np.matlib.repmat([3],6,5) # cf
    fts = np.matlib.repmat([0.80,0.65,0.50,0.35,0.20],6,1) # fts
    mle = np.matlib.repmat([4],6,5)
    trial = np.repeat([1, 2], 3)[:, None].repeat(5, axis=1)
    AMPO = dataExp[0:6,[5,6,2,7,8]]
    
    data['specimen'].extend([mus]*30)
    data['cf'].extend(cf.flatten())
    data['fts'].extend(fts.flatten())
    data['mle'].extend(mle.flatten())
    data['trial'].extend(trial.flatten())
    data['AMPO'].extend(AMPO.flatten())
    data['subset'].extend(['C']*30)

    # Panel D: 8mm MLE - effect of FTS
    cf = np.matlib.repmat([2],6,5) # cf
    fts = np.matlib.repmat([0.80,0.65,0.50,0.35,0.20],6,1) # fts
    mle = np.matlib.repmat([8],6,5)
    trial = np.repeat([1, 2], 3)[:, None].repeat(5, axis=1)
    AMPO = dataExp[6:12,[5,6,2,7,8]]
    
    data['specimen'].extend([mus]*30)
    data['cf'].extend(cf.flatten())
    data['fts'].extend(fts.flatten())
    data['mle'].extend(mle.flatten())
    data['trial'].extend(trial.flatten())
    data['AMPO'].extend(AMPO.flatten())
    data['subset'].extend(['D']*30)

df = pd.DataFrame(data)

#%% Create function to run mixedLM
def run_mixed_model(df, subset, predictor, xlabel=None):
    df_sub = df[df['subset'] == subset].copy()
    df_sub = df_sub.dropna(subset=['AMPO'])

    # Model fit
    formula = f"AMPO ~ {predictor} + I({predictor}**2)"
    model = mixedlm(formula, df_sub, groups=df_sub["specimen"])
    result = model.fit()

    #print(f"\n=== Subset {subset} ({predictor}) ===")
    #print(result.summary())
    
    # P-values
    pvals = result.pvalues
    p_lin = pvals[predictor]
    p_quad = pvals[f'I({predictor} ** 2)']
    
    if xlabel is not None:
        # Prediction
        x_vals = np.linspace(df_sub[predictor].min(), df_sub[predictor].max(), 100)
        x2_vals = x_vals**2
    
        fe = result.fe_params
        y_pred = (
            fe['Intercept']
            + fe[predictor]*x_vals
            + fe[f'I({predictor} ** 2)']*x2_vals
        )
        
        # Plot
        plt.figure(figsize=(6,4))
    
        for mus in df_sub['specimen'].unique():
            df_mus = df_sub[df_sub['specimen'] == mus]
            plt.scatter(df_mus[predictor], df_mus['AMPO'], label=mus, alpha=0.7)
    
        plt.plot(x_vals, y_pred, 'k-', linewidth=2, label='Mixed model fit')
    
        plt.xlabel(xlabel)
        plt.ylabel('AMPO')
        plt.title(f'Subset {subset}: AMPO vs {predictor}')
        plt.ylim(0, np.ceil(np.max(df_sub['AMPO'])/10)*10)
        plt.legend()
        plt.tight_layout()
        plt.show()

    return result, p_lin, p_quad

#%% Run linear mixed models
# Panel A & B: effect CF
result_A, pl_A, pq_A = run_mixed_model(df, 'A', 'cf') #, 'Cycle frequency [Hz]')
result_B, pl_B, pq_B = run_mixed_model(df, 'B', 'cf') #, 'Cycle frequency [Hz]')

# Panel C & D: effect FTS
result_C, pl_C, pq_C = run_mixed_model(df, 'C', 'fts') #, 'FTS')
result_D, pl_D, pq_D = run_mixed_model(df, 'D', 'fts') #, 'FTS')

# Store all p-values
pl = np.array([pl_A, pl_B, pl_C, pl_D]) # linear terms
pq = np.array([pq_A, pq_B, pq_C, pq_D]) # quadratic terms
p_all = np.array([pl, pq]) # all p-values

Compute differences in AMPO

"""
This script computes selected percentage differences in experimentally measured
AMPO for reporting in the manuscript.

The calculations summarise how AMPO changes with cycle frequency and FTS for
the 4 mm and 8 mm MTC length excursion conditions. Quadratic fits are used only
to estimate the approximate cycle frequency at which AMPO peaks for FTS = 0.5.
"""

#%% Load packages & set directories
import os, sys
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path

# Set directories
cwd = Path.cwd()
baseDir = cwd.parent.parent
dataDir = baseDir / 'data'
funcDir = baseDir / 'analysis' / 'functions'
sys.path.append(str(funcDir))

plt.close('all')

import stats

#%% Estimate optimum cycle frequency at 4 mm MTC length excursion
AMPO = np.empty((0,5))
for mus,clr in zip(['GMe1','GMe2','GMe3'],['k','r','g','b']):
    filepath = os.path.join(dataDir,mus,str(mus)+'_dataAMPO.xlsx')
    dataExp = pd.read_excel(filepath).to_numpy()[:,5:].astype(float)
    
    iCol = np.r_[1,2,3,4,5]-1 
    dataExp[dataExp<1] = np.nan
    AMPO = np.vstack((AMPO,dataExp[3:6,iCol]))

idx = np.isfinite(AMPO)
cf = np.tile([1,2,3,4,5],(9,1))
coef = np.polyfit(cf[idx].flatten(),AMPO[idx].flatten(),2)

fig, ax = plt.subplots()   
ax.plot(cf.T,AMPO.T)

cf = np.linspace(1,5,1000)
AMPOfit = np.polyval(coef,cf)
ax.plot(cf,AMPOfit,'--')

iMax = np.argmax(AMPOfit)
print(f"For FTS = 0.5 and MLE = 4 mm, optimum cycle frequency ≈ {cf[iMax]:.1f} Hz")

#%% Estimate optimum cycle frequency at 8 mm MTC length excursion
AMPO = np.empty((0,5))
for mus,clr in zip(['GMe1','GMe2','GMe3'],['k','r','g','b']):
    filepath = os.path.join(dataDir,mus,str(mus)+'_dataAMPO.xlsx')
    dataExp = pd.read_excel(filepath).to_numpy()[:,5:].astype(float)
    
    iCol = np.r_[1,2,3,4,5]-1 
    dataExp[dataExp<1] = np.nan
    AMPO = np.vstack((AMPO,dataExp[9:12,iCol]))

idx = np.isfinite(AMPO)
cf = np.tile([1.0,1.5,2.0,2.5,3.0],(9,1))
coef = np.polyfit(cf[idx].flatten(),AMPO[idx].flatten(),2)

fig, ax = plt.subplots()   
ax.plot(cf.T,AMPO.T)

cf = np.linspace(1,3,1000)
AMPOfit = np.polyval(coef,cf)
ax.plot(cf,AMPOfit,'--')

iMax = np.argmax(AMPOfit)
print(f"For FTS = 0.5 and MLE = 8 mm, optimum cycle frequency ≈ {cf[iMax]:.1f} Hz")

#%% Compute FTS effects at 4 mm MTC length excursion
print()
print('MTC length excursion = 4 mm')

AMPO = np.empty((0,13))
for mus,clr in zip(['GMe1','GMe2','GMe3'],['k','r','g','b']):
    filepath = os.path.join(dataDir,mus,str(mus)+'_dataAMPO.xlsx')
    dataExp = pd.read_excel(filepath).to_numpy()[:,5:].astype(float)
    
    dataExp[dataExp<1] = np.nan
    AMPO = np.vstack((AMPO,dataExp[3:6,:]))

# FTS 0.20 -> 0.50 - 4 mm @ 3Hz
pDiff = np.nanmean(stats.pdiff(AMPO[:,2],AMPO[:,8]))
print(f"FTS 0.20 -> 0.50 - 4mm@3Hz: {pDiff:.1f} %")

# FTS 0.50 -> 0.80 - 4 mm @ 3Hz
pDiff = np.nanmean(stats.pdiff(AMPO[:,5],AMPO[:,2]))
print(f"FTS 0.50 -> 0.80 - 4mm@3Hz: {pDiff:.1f} %")

# FTS 0.65 -> 0.80 - 4 mm @ 3Hz
pDiff = np.nanmean(stats.pdiff(AMPO[:,5],AMPO[:,6]))
print(f"FTS 0.65 -> 0.80 - 4mm@3Hz: {pDiff:.1f} %")

# FTS 0.50 -> 0.80 - 4 mm @ 5Hz
pDiff = np.nanmean(stats.pdiff(AMPO[:,9],AMPO[:,4]))
print(f"FTS 0.50 -> 0.80 - 4mm@5Hz: {pDiff:.1f} %")

#%% Compute FTS effects at 8 mm MTC length excursion
print()
print('MTC length excursion = 8 mm')
AMPO = np.empty((0,13))
for mus,clr in zip(['GMe1','GMe2','GMe3'],['k','r','g','b']):
    filepath = os.path.join(dataDir,mus,str(mus)+'_dataAMPO.xlsx')
    dataExp = pd.read_excel(filepath).to_numpy()[:,5:].astype(float)
    
    dataExp[dataExp<1] = np.nan
    AMPO = np.vstack((AMPO,dataExp[9:12,:]))

# FTS 0.20 -> 0.50 - 8 mm @ 2Hz
pDiff = np.nanmean(stats.pdiff(AMPO[:,2],AMPO[:,8]))
print(f"FTS 0.20 -> 0.50 - 8mm@2Hz: {pDiff:.1f} %")

# FTS 0.50 -> 0.80 - 8 mm @ 2Hz
pDiff = np.nanmean(stats.pdiff(AMPO[:,5],AMPO[:,2]))
print(f"FTS 0.50 -> 0.80 - 8mm@2Hz: {pDiff:.1f} %")

# FTS 0.65 -> 0.80 - 8 mm @ 2Hz
pDiff = np.nanmean(stats.pdiff(AMPO[:,5],AMPO[:,6]))
print(f"FTS 0.65 -> 0.80 - 8mm@2Hz: {pDiff:.1f} %")

# FTS 0.50 -> 0.80 - 4 mm @ 3Hz
pDiff = np.nanmean(stats.pdiff(AMPO[:,9],AMPO[:,4]))
print(f"FTS 0.65 -> 0.80 - 8mm@3Hz: {pDiff:.1f} %")