Enhancing spontaneous recovery after stroke: a randomized controlled trial.

1. % All-in-one pipeline that integrates signal processing, data analysis, data formatting, and S-R curve fitting.
3. % A Matlab sigmoid function is fitted which returns 4 parameters [minimum y value, maximum y value, S50, slope parameter k]
4. % This code has adjustable constraints to improve the fit of the curve (e.g., y values cannot be negative).
5. % There are also adjustable inital parameters (search 'initialparams') to give the algorithm somewhere to start.
7. %% LOADING THE FILES
8. % Currently this code runs one participants data at one time point.
9. % All variables need to be cleared in between running data, the code does this in the very first line. If there's just one section you want to run
10. % again (e.g., SR curve section with altered parameters) clear all and run again from the start.
11. %
12. % The code requires a signal file in .mat format. All of the matlab files are on the ESPRESSo research drive and are ready to be inputted to this code and processed.
13. % The name of the file gets entered in the file prefix section in General Info.
14. % The required input is the file prefix. Please enter the ESPRESSo alphanumeric code
15. % Please make sure to exclude the three zeros at the end.
16. % E.g., '#########_SR_L_XX##ESP'
17. %
18. %% MODIFIABLE ASPECTS OF THE CODE
19. % Enter your initials in the assessor variable in the General Info section (search 'assessor').
20. %
21. % The second argument asks for the side that was analysed (specify as either 'P/NP').
22. %
23. % If the SR fit is not great the first time around, change the initial parameters (search 'initialparams' within the script to find it) to give a more accurate prediction of S50 and/or k.
24. %
25. % The title of the S-R curve will need to change between P/NP (search 'title')
26. %
27. %% SAVING THE OUTPUT FILES
28. % The output file and SR curve image will automatically be saved in an Output File folder within the folder that you are using the data from.
29. % The folder will have the initials of the assessor at the start of it.
30. % Enter your initials in the assessor variable in the General Info section.
31. %
32. % When you press the Run button you will have the SR curve in a separate figure window and an Output structure and SlopeAtS50 (in mV/1%MSO) variable in the
33. % workplace.
35. % SR CURVE FIGURE
36. % The SR curve will pop up in a separate figure. It will the plots for the 4 muscles in different colours and a legend to label the muscles.
37. % It will have shaded regions at the bottom of the plot showing the background RMS lower and upper 95% confidence interval. The shaded regions are colour coded to match the muscle.
38. % The x-axis is fixed from 0 to 100 and the y-axis is auto scaled.
39. %
40. % OUTPUT STRUCTURE IN WORKSPACE
41. % Inside the output structure are all the variables calculated including the slope variables, I've just put a copy outside the structure
42. % too for easy access. The output structure is what get saved for each participant.
43. % The conventional reporting for the outputs that do not have the muscle names explicitly displayed is:
44. % 1 = FDI
45. % 2 = ADM
46. % 3 = ECR
47. % 4 = FCR
48. % The output structure contains:
49. % FDI, ADM, ECR, FCR structures - all the information pertaining to just the specfic muscle
50. % Data table
51. % State and stimulus intensity (linked i.e., state 1 = 30% MSO, state 2 = 40% MSO etc.)
52. % BGamp = background amplitude calculated in the prestimulus window
53. % defined in the script (-55 to -5 ms unless changed)
54. % MEPamp = average amplitude for each state/intensity. Note the state 0
55. % value is the average BGamplitude value. This is to help with curve
56. % fitting to anchor the curve to a value at x = 0.
57. % StD = standard deviation of the MEP amplitudes
58. % RMS = background RMS calculated in the prestimulus window defined in
59. % the script (-55 to -5 ms unless changed)
60. % Norm2Max = MEP amplitude values normalised to the maximum MEP
61. % amplitude value. Not being used for analysis currently but retained
62. % it in the output just in case
63. % SI = array of the stimulus intensities used for SR curve (i.e.,
64. % including 0% MSO)
65. % MEPvalues table = actual MEP amplitude values used for SR curve
66. % SigX = array of x values plotted for the actual sigmoid function
67. % SigY = array of y values plotted for the sigmoid function
68. % S50 = S50 value (i.e. x value when y value is half it's maximum)
69. % k = slope parameter value k. Used to calculate slope at S50 with the
70. % formula Max/4*k
71. % Slope_at_S50 = slope of the curve at S50 in mV per 10% MSO
72. % The following tables in the Output are contain the same information as what is described in the individual muscle structures. It is just displayed in a table showing all 4 muscles at once.
73. % MEPamplitude table
74. % StD table
75. % BGamp table
76. % RMS table
77. % RMSlowerCI = used as the lower boundary of the shaded area of RMS on
78. % the SR curve
79. % RMSupperCI = used as the upper boundary of the shaded area of RMS on
80. % the SR curve
81. % AverageRMS = one average background RMS value for each muscle across
82. % all trials
83. % RMS_thresholds = the lowest background rms threshold (starting at 10
84. % microvolts) for each muscle to ensure at least 5 frames at each intensity
85. % k = slope parameter k value for the 4 muscles
86. % S50 = S50 value for all four muscles (i.e. x value when y value is half it's maximum)
87. % SlopeAtS50_mV1MSO = the slope of the SR curve at S50 for all 4
88. % muscles in mV per 1%MSO
89. % SlopeAtS50_mV10MSO = the slope of the SR curve at S50 for all 4
90. % muscles in mV per 10%MSO
91. % GoodnessofFit = Goodness of fit of the curve to the data points using
92. % total sums squared and r squared values
93. % FrameNumber table = showing the number of frames used for each muscle
94. % at each intensity after excluding based on background RMS. This will
95. % always be at least 5.
96. % TotalTrials table = showing the total number of frames for each
97. % muscles/intensity combination before any cleaning/trimming of the
98. % data
99. % MEPsum table = the sum of the MEP values used to create the S-R curve
100. % for each muscle (i.e., the sum of the average MEP amplitude at each
101. % stimulus intensity)
103. % Code created by Maxine Shanks and Gene Tang
104. % Last update: Maxine Shanks 29-08-2024 with the update of checking that the
105. % estimated RMT matches the actualy datapoints and adjusting accordingly if not
107. %% General Info
109. clear all;
110. close all;
111. assessor = 'MS' % Enter your initials
112. n_channels = 4; % set number of channels/muscles recorded from
113. fileprefix = '180324_SR_L_xxxxESP' % enter the name of the file to analyse (excluding the part of the filename that is in line 31 e.g., '000.mat')
114. side = 'NP' % enter the label for the side to be analysed (e.g., L/R, D/ND, P/NP)
115. Target_muscles_labels = {'FDI' 'ADM' 'ECR' 'FCR'}; % enter the names for the 4 muscle labels
116. channelNames = {[side Target_muscles_labels{1}],[side Target_muscles_labels{2}],[side Target_muscles_labels{3}],[side Target_muscles_labels{4}]}; % the channel names will then become the side+channel name e.g., PFDI for paretic FDI
117. fs = 2000; % sampling rate of the signal
119. %% Loading data and assigning time
120. r_data = load([fileprefix '000.mat']); % load in rest data using the fileprefix set above in General Info + the suffix i.e., '000.mat'
121. fnames = fieldnames(r_data); % extract relevant data from imported .cfs file
122. r_data = r_data.(fnames{1});
124. % Create table to store relevant information in
125. restData = struct2table(r_data.frameinfo); % put frameinfo into a table
126. restData.points = []; % remove unwanted points column
127. restData.sweeps = []; % remove unwanted sweeps column
129. % Create list of times from -0.5 to 0.5 to correspond with the samples
130. times = [];
131. for i = 1:length(r_data.values)
132. times(i) = i *(1/fs)-0.5; % convert times to -0.5 to 0.5 pre/post trigger
133. end
135. % Set background analysis window
136. baseline_time = [-0.055, -0.005]; % use time point -50ms to -5ms to calculate background rms
137. baseline_s = find(times >= baseline_time(1), 1):find(times <= baseline_time(2), 1, 'last'); % find times that are within the background window
139. % Set p2p time
140. MEP_time = [0.015, 0.06]; % setting MEP window to 15-60 ms
141. MEP_s = find(times >= MEP_time(1), 1):find(times <= MEP_time(2), 1, 'last'); % find times that are within the p2p window
143. % Create temporary data to work on
144. temp_signal_data = r_data;
145. temp_data = restData;
147. % compute RMS/RMS amplitude/MEP amplitude/MEP latency
148. for c = 1:n_channels % iterate through the number of channels there are (e.g., 4)
150. for f = 1:size(temp_data.number) % iterate through the number of frames (e.g., 90)
152. % Background RMS
153. temp_data.rms(f,:) = rms(r_data.values(baseline_s(1,1):baseline_s(1,end),c,f)); % calculate rms in the background window set above (e.g., -50 to -5 ms)
155. % Background amplitude
156. temp_data.BGamp(f,:) = peak2peak(r_data.values(baseline_s(1,1):baseline_s(1,end),c,f)); % caluculate the peak2peak amplitude in the background window
158. % MEP amplitude
159. temp_data.MEPamp(f,:) = peak2peak(r_data.values(MEP_s(1,1):MEP_s(1,end),c,f)); % calculate the peak2peak amplitude in the p2p window set above (e.g., 0.01 to 0.06)
161. end
163. % Update variable names in data table
164. temp_data.Properties.VariableNames{'rms'} = [channelNames{c} '_rms']; % creates a rms column with the data calculated above
165. temp_data.Properties.VariableNames{'BGamp'} = [channelNames{c} '_BGamp']; % creates a Bgamp column with the data calculated above
166. temp_data.Properties.VariableNames{'MEPamp'} = [channelNames{c} '_MEPamp']; % creates a MEPamp column with the data calculated above
168. end
170. % Store temporary data back into relevant variables
171. raw_data = temp_data;
172. raw_data.start = [];
173. raw_data = sortrows(raw_data,'state','ascend'); % sort the 90 frames by state in ascending order (e.g., 30% - 100% MSO)
174. rest_signal_data = temp_signal_data;
176. n_frames = size(raw_data,1);
178. clear baseline_s f i MEP_s rest_signal_data restData temp_data temp_signal_data % clears all the unnecessary variables from the workspace
179. %% Cleaning the data
180. %this section takes BGamp MEPamp and rms of each muscle and assigns them
181. %into respective variable.
183. % Splits the background rms, background amplitude and MEP amplitude data into the individual muscles
184. State_index = raw_data(:,"state");
185. FDI_temp= raw_data(:,[{'state'}, {[side 'FDI_BGamp']}, {[side 'FDI_MEPamp']}, {[side 'FDI_rms']}]);
186. ADM_temp = raw_data(:,[{'state'}, {[side 'ADM_BGamp']}, {[side 'ADM_MEPamp']}, {[side 'ADM_rms']}]);
187. ECR_temp = raw_data(:,[{'state'}, {[side 'ECR_BGamp']}, {[side 'ECR_MEPamp']}, {[side 'ECR_rms']}]);
188. FCR_temp = raw_data(:,[{'state'}, {[side 'FCR_BGamp']}, {[side 'FCR_MEPamp']}, {[side 'FCR_rms']}]);
190. var_names = {'FDI_temp', 'ADM_temp', 'ECR_temp', 'FCR_temp'}; % sets the names of the 4 muscles
191. outliers_flags = []; % creates an array to store outlier flags. %matrix contains 4 columns each for each muscle. 1 represents the row that an outlier has been removed
192. var_names = {'FDI_temp', 'ADM_temp', 'ECR_temp', 'FCR_temp'};
193. muscle_names = {'FDI', 'ADM', 'ECR', 'FCR'};
194. remaining_trials_percentage = 0.3;
195. threshold = 0.5;
196. rms_threshold_default = 0.010; % Initial RMS threshold in mV
197. rms_thresholds = rms_threshold_default * ones(1, numel(var_names)); % Initialize RMS thresholds for each muscle
198. frames_per_intensity = zeros(9, numel(var_names)); % Initialize the array to store frame counts
199. total_trials = ones(9, numel(var_names));
200. outliers_flags = false(size(eval(var_names{1}), 1), numel(var_names)); % Initialize the array to document outliers
201. original_data = struct(); % Create a structure to store the original data
203. for vx = 1:numel(var_names) % cycle through the muscles
204. original_data.(var_names{vx}) = eval(var_names{vx}); % grab the data for the muscle
205. end
207. % Loop through all states and make sure there are at least the minimum number of frames remaining after the exclusion based on background RMS.
208. % The RMS criterion increases by 0.005 mV each iteration if there's not enough frames across all intensities. This process works for each muscle
209. % separately
210. while any(frames_per_intensity < threshold*total_trials, 'all')%remaining_trials_percentage < 0.5)
211. for vx = 1:numel(var_names) % For each muscle
212. var_value = original_data.(var_names{vx}); % Restore the original _temp data;
213. rms_threshold = rms_thresholds(vx);
215. for RMS_col = 4 % For the RMS column
216. outlier_index = (var_value{:,RMS_col}) > rms_threshold;
217. outliers_flags(:, vx) = outliers_flags(:, vx) | outlier_index; % Store the outlier index into a column corresponding to the muscle
219. var_value{outlier_index, 2:end} = NaN; % Turn all the rows with outliers into NaNs
220. eval([var_names{vx} ' = var_value']);
221. end
222. end
224. for vx = 1:numel(var_names) % For each muscle
225. var_value = eval(var_names{vx}); % Get the muscle's _temp table
226. increase_threshold = true; % Flag to indicate whether threshold needs to be increased
228. for state = 1:9 % For each state
229. valid_rows = sum(~isnan(var_value{var_value.state == state, 2}), 2); % Count valid rows for the state and store as 1's. Non valid frames are 0.
230. total_frames = sum(valid_rows); % total number of valid frames
231. frames_per_intensity(state,vx) = total_frames;
233. total_trials(state, vx) = sum(var_value.state == state);
235. min_trials_required = ceil(total_trials(state,vx)*0.5);
237. remaining_trials_percentage = total_frames/total_trials(state,vx);
239. if remaining_trials_percentage < 0.5 && increase_threshold
240. rms_thresholds(vx) = rms_thresholds(vx) + 0.005;
241. increase_threshold = false;
242. end
243. end
244. end
245. end
247. avg_rms_values = zeros(1, numel(var_names)); % Initialize the array to store average RMS values
249. for vx = 1:numel(var_names) % For each muscle
250. var_value = eval(var_names{vx}); % Get the muscle's _temp table
251. rms_column = var_value{:, 4}; % Get the RMS values from column 4
253. % Calculate the average RMS value excluding NaNs
254. avg_rms_values(vx) = mean(rms_column, 'omitnan');
255. end
257. outliers_flags = array2table(outliers_flags, 'VariableNames', {[side 'FDI'], [side 'ADM'], [side 'ECR'], [side 'FCR']}); % create a table for outliers_flags
259. data = horzcat(State_index, FDI_temp(:,2:end), ADM_temp(:,2:end), ECR_temp(:,2:end), FCR_temp(:,2:end)); % concatenate all tables back into 1 file called rest data
261. clear outlier_index outliers_flags RMS_col State_index vx var_value var_names increase_threshold rms_column rms_threshold_default rms_threhsold total_frames valid_rows state min_frames_per_intensity
262. %% Calculate percentage of frames retained
264. frames_retained = sum(sum(frames_per_intensity)); %add all elements of the array
265. percent_retained = (frames_retained/(n_frames*4))*100 %calculte percentage of frames retained out of 360 total frames
266. %% Calculate averages of BGamp/MEPamp/RMS
268. BGamp = table; %create table for each values
269. MEPamp = table;
270. RMS = table;
271. StD = table;
272. SI = [30;40;50;60;65;70;80;90;100]; % set the intensities
273. BGamp.Intensities = SI; % assign the above intensities to the following tables
274. MEPamp.Intensities = SI;
275. RMS.Intensities = SI;
276. StD.Intensities = SI;
277. n_states = 9; % number of states
278. n_muscles = 4; % number of muscles
280. for c = 1:n_channels % cycle through each of the muscles/channels
281. var = [channelNames{c} '_BGamp']; % the variable is set as the BGamp column
282. tidx = find(string(data.Properties.VariableNames) == var); % find the BGamp column
283. for s = 1:n_states % for each of the states
284. for f = 1:size(data(data.state==s,:),1) % for each of the frames
285. state_data = data(data.state==s,:);
286. end
287. BGav(s) = mean(table2array(data(data.state==s, tidx)), 'omitnan'); % calculate the mean of the BGamplitude for each state, omitting any NaN values
288. MEPav(s) = mean(table2array(data(data.state==s, tidx+1)), 'omitnan'); % calculate the mean of the MEPamplitude for each state, omitting any NaN values
289. RMSav(s) = mean(table2array(data(data.state==s, tidx+2)), 'omitnan'); % calculate the mean of the RMS for each state, omitting any NaN values
290. Stdev(s) = std(table2array(data(data.state==s, tidx+1)), 'omitnan'); % calculate the standard deviation of the MEPamplitude values for each state, omitting any NaN values
292. end
294. BGamp.([channelNames{c}]) = BGav([1:9]).'; % create a table for the background amplitude with the average value for each muscle at each intensity
295. MEPamp.([channelNames{c}]) = MEPav([1:9]).'; % create a table for the MEP amplitude with the average value for each muscle at each intensity
296. RMS.([channelNames{c}]) = RMSav([1:9]).'; % create a table for the RMS with the average value for each muscle at each intensity
297. StD.([channelNames{c}]) = Stdev([1:9]).'; % create a table for the St Dev with the average value for each muscle at each intensity
299. end
301. clear BGav c f MEP_time MEPav n_channels n_muscles n_states RMSav s state_data Stdev tidx var
302. %% Reformat the data for S-R curve plotting
303. %this section takes BGamp MEPamp and RMS values from the previous section
304. %and assign them into their respective muscles (table format)
306. FDI = table; ADM = table; ECR = table; FCR = table; % create a table for each muscle
307. Target_muscles_table = {FDI ADM ECR FCR};
309. State = [1;2;3;4;5;6;7;8;9]; % set the state names
310. SI = [30;40;50;60;65;70;80;90;100]; % set the SI column
312. for mx = 1:length(Target_muscles_labels) % for each of the 4 muscles
313. Target_muscles_table{mx}.State = State; % assigning state array to the state column
314. Target_muscles_table{mx}.SI = SI; % assigning intensity to the SI column
315. Target_muscles_table{mx}.BGamp = table2array(BGamp(:,[side,Target_muscles_labels{mx}])); % assinging BGamp to BGamp column
316. Target_muscles_table{mx}.MEPamp = table2array(MEPamp(:,[side,Target_muscles_labels{mx}])); % assinging MEPamp to MEPamp column
317. Target_muscles_table{mx}.StD = table2array(StD(:,[side,Target_muscles_labels{mx}])); % assigning StDev to StDev column
318. Target_muscles_table{mx}.RMS = table2array(RMS(:,[side,Target_muscles_labels{mx}])); % assinging RMS to RMS column
319. norm2max = [];
320. for sx = 1:9
321. norm2max(sx,1) = table2array(Target_muscles_table{mx}(sx,4))/(max(Target_muscles_table{mx}{:,4})); % creating a normalised to max value (currently not using)
322. end
323. Target_muscles_table{mx}.Norm2Max = norm2max; % adding normalised values to the table (currently not using)
324. end
326. clear sx State norm2max mx
327. %% Plot the points and add function stats to the table
328. %This section plots x = SI, y = MEPamp and fit a sigmoid function.
330. colors = {'#A2142F', '#0072BD', '#77AC30', '#D95319'}; % specify colour of the plot
332. Target_muscles_labels = {'FDI' 'ADM' 'ECR' 'FCR'};
334. % Temporarily storing S50 and slope of all 4 muscles
335. SIarray = zeros(4,10); % array to store the SI data points
336. MEPampArray = zeros(4,10); % array to store the MEPamp data points
337. SigXarray = zeros(4,101); % array to store the sigmoid x values
338. SigYarray = zeros(4,101); % array to store the sigmoid y values
339. S50 = []; % array to store S50
340. k = []; % array to store the slope parameter
341. dy = {}; % cell to contain the derivative of the sigmoid functions
342. Max = []; % storing max slope from the derivative
343. RMSmean = zeros(10,4); % array to store the RMS mean
344. RMSupperCI = zeros(10,4); % array to store the RMS upper confidence interval
345. RMSlowerCI = zeros(10,4); % array to store the RMS lower confidence interval
346. GoodnessofFit = zeros(1,4); % array to store the goodness of fit of the SR curve
347. estimated_threshold = zeros(1,4); % array to store the estimated RMT from the sigmoid curve
348. adj_estimated_threshold = zeros(1,4); % array to store the adjusted estimated RMT from the sigmoid curve
349. new_adj_estimated_threshold = zeros(1,4); % array to store the new adjusted estimated RMT from the actual datapoints if the eRMT from the sigmoid doesn't make sense
351. fdi_bgampstd = std(table2array(BGamp(:,2))); % finding the standard deviation of the background amplitude for FDI
352. adm_bgampstd = std(table2array(BGamp(:,3))); % finding the standard deviation of the background amplitude for ADM
353. ecr_bgampstd = std(table2array(BGamp(:,4))); % finding the standard deviation of the background amplitude for ECR
354. fcr_bgampstd = std(table2array(BGamp(:,5))); % finding the standard deviation of the background amplitude for FCR
355. newRowStD = [0,fdi_bgampstd,adm_bgampstd,ecr_bgampstd,fcr_bgampstd]; % creating a new array with the 4 st deviation values
356. newRowStD = array2table(newRowStD); % convert the new array to a table
357. newRowStD.Properties.VariableNames = {'Intensities',channelNames{1,1},channelNames{1,2},channelNames{1,3},channelNames{1,4}}; % Naming the variables in the table
358. StD = [newRowStD; StD]; % adding the new standard deviation table to the existing (i.e., adding in a 0% MSO row)
360. for ii = 1:numel(Target_muscles_labels) %loop through each muscle
361. values = Target_muscles_table{ii};
363. bgampav = mean(table2array(values(:,3))); % finding the mean background amplitude value
364. newRowAmp = [0,0,0,bgampav,table2array(StD(1,ii+1)),0,0]; % creating a new array with the mean background amplitude
365. newRowAmp = array2table(newRowAmp); % convert the new array to a table
366. newRowAmp.Properties.VariableNames = {'State', 'SI', 'BGamp', 'MEPamp', 'StD', 'RMS', 'Norm2Max'}; % Naming the variables in the table
367. values = [newRowAmp; values]; % adding the new background amplitude row to the existing table (i.e., 0% MSO row)
369. % Calculating the 95% CI for RMS from the state 1 (i.e., 30% MSO)
370. % for the values to use at 0% MSO on the graph
371. muscle_name = Target_muscles_labels{ii}; % Get muscle name
372. muscle_temp_table_name = [muscle_name '_temp']; % Construct the table name
373. muscle_temp_table = eval(muscle_temp_table_name); % access the specific muscle's _temp table
374. alpha = 0.05; % significance level for 95% confidence level
375. valid_rows = sum(~isnan(muscle_temp_table{muscle_temp_table.state == 1, 4}), 2); % Get the valid rows of the RMS column for state 1 (30% MSO)
376. n = sum(valid_rows); % number of valid rows
377. mean_RMS_state1 = mean(muscle_temp_table{muscle_temp_table.state == 1,4}, 'omitnan'); % find the mean, omitting NaN values
378. std_RMS_state1 = std(muscle_temp_table{muscle_temp_table.state == 1,4}, 'omitnan'); % find the standard deviation, excluding NaN values
379. t_score = tinv(1-alpha/2, n-1);
380. margin_of_error = t_score*(std_RMS_state1/sqrt(n));
381. CI_lower_BGamp = mean_RMS_state1 - margin_of_error;
382. CI_upper_BGamp = mean_RMS_state1 + margin_of_error;
383. RMSmean(1,ii) = mean_RMS_state1; % storing the mean RMS value
384. RMSlowerCI(1,ii) = CI_lower_BGamp; % storing the RMS lower CI
385. RMSupperCI(1,ii) = CI_upper_BGamp; % storing the RMS upper CI
387. % Calculating the 95% Confidence Interval for RMS for the states
388. % (i.e., 30 - 100% MSO). Same process as above but cycling through
389. % the states
390. for s = 1:numel(unique(muscle_temp_table.state))
391. rows_of_interest = muscle_temp_table.state == s;
392. valid_rows = sum(~isnan(muscle_temp_table{rows_of_interest, 4}), 2);
393. n = sum(valid_rows);
394. mean_rms_for_state = mean(muscle_temp_table{rows_of_interest, 4}, 'omitnan');
395. std_rms_for_state = std(muscle_temp_table{rows_of_interest, 4}, 'omitnan');
396. t_score = tinv(1-alpha/2, n-1);
397. margin_of_error = t_score*(std_rms_for_state/sqrt(n));
398. CI_lower = mean_rms_for_state - margin_of_error;
399. CI_upper = mean_rms_for_state + margin_of_error;
400. RMSmean(s+1,ii) = mean_rms_for_state;
401. RMSlowerCI(s+1,ii) = CI_lower;
402. RMSupperCI(s+1,ii) = CI_upper;
403. end
405. % Assigning the x and y values
406. x = table2array(values(:,"SI")); %assiging stim intensity as X
407. y = table2array(values(:,"MEPamp")); %assigning MEPamp as y
409. % Define sigmoid function
410. sigmoid = @(params, x) params(1) + (params(2) - params(1)) ./ (1 + exp(-(x - params(3)) / params(4))); % sigmoid curve equation
412. % Initial parameter guess [min y value, max y value, S50, k] A good
413. % starting value for S50 is 65 and k of 5
414. initialParams = [min(y),max(y),65,5]; % Having the 0 for min(y) is ensure min(y) is non-negative. These initial parameters can be changed to refine the fit.
416. % Fit the sigmoid function to the data using lsqcurvefit
417. fitParams = lsqcurvefit(sigmoid, initialParams, x, y);
419. % Additional fmincon optimization for parameter constraints
420. lb = [0,0,20, 1]; % lower bounds for parameters
421. ub = [min(y),15,90,20]; % upper bounds for parameters
422. options = optimoptions('fmincon', 'Display', 'off'); % Set optimization options
423. fitParams_constrained = fmincon(@(params)sum((y - sigmoid(params, x)).^2), fitParams, [], [], [], [], lb, ub, [], options);
425. % Generate curve using the fitted parameters
426. curveX = 0:100;
427. curveY_constrained = sigmoid(fitParams_constrained, curveX);
429. % Plot data and fitted curve for the current muscle
430. plot(x, y, 'Color', colors{1,ii}, 'Marker', '.', 'LineStyle', 'none', 'MarkerSize', 10,'HandleVisibility','off'); % Plot data for the current muscle
431. hold on;
432. plot(curveX, curveY_constrained, 'Color', colors{1,ii}, 'DisplayName', Target_muscles_labels{ii}, 'LineWidth', 2);
433. hold off;
435. Target_muscles_table{ii} = values;
437. %Store S50 and slope for each muscle in a 1x4 array
438. SIarray(ii,:) = x;
439. MEPampArray(ii,:) = y;
440. SigXarray(ii,:) = curveX;
441. SigYarray(ii,:) = curveY_constrained;
442. S50(1,ii) = fitParams_constrained(1,3); % S50 the 3rd value from the params output
443. k(1,ii) = (fitParams_constrained(1,4)); %divide by 10 to adjust for magnitude scaling %no longer used
444. SlopeAtS50(1,ii) = fitParams_constrained(1,2)/(4*fitParams_constrained(1,4)); % equation for getting slope parameter (k) into slope at S50. Equation = Max/4*k
446. %plot error bars
447. % hold on
448. % errorbar(x,y,table2array(StD(:,ii+1)));
449. % eb = get(gca, 'Children');
450. % set(eb(1), 'Color', colors{ii}, 'DisplayName', Target_muscles_labels{ii},'Linewidth', 1,'HandleVisibility','off','Linestyle','none');
452. % Convert hexidecimal colour to RGB triplet
453. hex_color = colors{ii};
454. color_rgb = sscanf(hex_color(2:end), '%2x%2x%2x', [1 3])/255;
455. % Plot the shaded region for the 95% CI for RMS
456. hold on;
457. fill([x; flipud(x)], [RMSlowerCI(:,ii); flipud(RMSupperCI(:,ii))], color_rgb, 'FaceAlpha', 0.5, 'EdgeColor', 'none','DisplayName', "Background RMS");
459. % Calculating TSS
460. % Calculate the mean of the y values (actual data points)
461. meanY = mean(y);
462. % Calculate the total sum of squares
463. tss = sum((y-meanY).^2);
464. % Matching the SigYarray values to the actual data points
465. indices = [1, 30, 40, 50, 60, 65, 70, 80, 90, 100];
466. fittedValues = SigYarray(ii,indices);
467. % Calculate the residuals
468. residuals = y.' - fittedValues;
469. % Calculate the sum of squares of the residuals (sum of squared differences between fitted curve and data points)
470. ssr = sum((residuals).^2);
471. % Calculate the coefficient of determination (R-squared value)
472. rSquared = 1 - (ssr/tss);
473. GoodnessofFit(1,ii) = rSquared;
475. % Identify the x value when the y value first goes above 0.05 for 3 data points
476. y_threshold = 0.05; % set the MEP amplitude value in mV for the threshold
477. consecutive_count = 0;
478. for j = 1:numel(SigXarray(1,:)) % cycle through the stimulus intensities in the Sigmoid curve fit
479. if SigYarray(ii,j) > y_threshold % if the MEP amplitude value from the sigmoid curve data points is above the threshold add one value to the consecutive count
480. consecutive_count = consecutive_count + 1;
481. if consecutive_count == 3 % once the consecutive count = 3
482. estimated_threshold(1,ii) = SigXarray(ii,j-(consecutive_count-1)); % go back to the first x value when y > 0.05 and put that in the estimated_threshold variable
483. break; % Stop searching when you find the x value
484. end
485. else
486. consecutive_count = 0; % Reset consecutive count if y is below the threshold
487. end
488. end
490. % Check if the threshold was not reached in the loop
491. if consecutive_count < 3
492. adj_estimated_threshold(1,ii) = 101; % Set to 101 if threshold is not reached
493. end
495. % Identify the x value when the y value first goes above 0.05 for 3 data points
496. adj_y_threshold = mean(SigYarray(ii,1:5)) + 0.05; % set the MEP amplitude value in mV for the threshold
497. consecutive_count2 = 0;
498. for j = 1:numel(SigXarray(1,:)) % cycle through the stimulus intensities in the Sigmoid curve fit
499. if SigYarray(ii,j) > adj_y_threshold % if the MEP amplitude value from the sigmoid curve data points is above the threshold add one value to the consecutive count
500. consecutive_count2 = consecutive_count2 + 1;
501. if consecutive_count2 == 3 % once the consecutive count = 3
502. adj_estimated_threshold(1,ii) = SigXarray(ii,j-(consecutive_count2-1)); % go back to the first x value when y > 0.05 and put that in the estimated_threshold variable
503. break; % Stop searching when you find the x value
504. end
505. else
506. consecutive_count2 = 0; % Reset consecutive count if y is below the threshold
507. end
508. end
510. % Check if the threshold was not reached in the loop
511. if consecutive_count2 < 3
512. adj_estimated_threshold(1,ii) = 101; % Set to 101 if threshold is not reached
513. end
515. % Find the closest x value upwards
516. closest_x_idx = find(x >= adj_estimated_threshold(1,ii),1);
517. closest_x = x(closest_x_idx);
519. % Get the corresponding y value
520. corresponding_y = y(closest_x_idx);
522. % Compare y value with adj_y_threshold
523. if corresponding_y < adj_y_threshold
524. % Move upwards until a y value is found that meets the criteria
525. for t = closest_x_idx+1:length(x)
526. if y(t) >= adj_y_threshold
527. closest_x_idx = t;
528. closest_x = x(t); % update closest x value to the new one
529. break; % stop the loop once the condition is met
530. end
531. end
532. % Take the value one lower than closest_x
533. new_eRMT = x(closest_x_idx - 1);
534. new_adj_estimated_threshold(1,ii) = new_eRMT;
535. else
536. new_adj_estimated_threshold(1,ii) = adj_estimated_threshold(1,ii);
537. end
539. % Ensure the minimum value of new_eRMT is not < 20% MSO
540. new_adj_estimated_threshold(1,ii) = max(new_adj_estimated_threshold(1,ii), 20);
542. end
544. % Creating the RMS CI values from S-R curve into a table
545. RMSlowerCI = array2table(RMSlowerCI); % convert the new array to a table
546. Intensities = [0;30;40;50;60;65;70;80;90;100];
547. RMSlowerCI = addvars(RMSlowerCI, Intensities, 'Before', 1);
548. RMSlowerCI.Properties.VariableNames = {'Intensities',channelNames{1,1},channelNames{1,2},channelNames{1,3},channelNames{1,4}}; % Naming the variables in the table
549. RMSupperCI = array2table(RMSupperCI);
550. RMSupperCI = addvars(RMSupperCI, Intensities, 'Before', 1);
551. RMSupperCI.Properties.VariableNames = {'Intensities',channelNames{1,1},channelNames{1,2},channelNames{1,3},channelNames{1,4}}; % Naming the variables in the table
552. RMSmean = array2table(RMSmean);
553. RMSmean = addvars(RMSmean, Intensities, 'Before', 1);
555. frames_per_intensity = array2table(frames_per_intensity);
556. frames_per_intensity = addvars(frames_per_intensity, SI, 'Before', 1);
557. frames_per_intensity.Properties.VariableNames = {'Intensities',channelNames{1,1},channelNames{1,2},channelNames{1,3},channelNames{1,4}}; % Naming the variables in the table
558. total_trials = array2table(total_trials);
559. total_trials = addvars(total_trials, SI, 'Before', 1);
560. total_trials.Properties.VariableNames = {'Intensities',channelNames{1,1},channelNames{1,2},channelNames{1,3},channelNames{1,4}};
562. % Add labels and title for the plot
563. xlim([0 100]); %set x-axis limit from 0 - 100
564. ylim([0,6]); %auto set the y-axis limit
565. xticks(0:10:100);
566. yticks(0:1:6);
567. xlabel('Stimulus intensity (%MSO)');
568. ylabel('Evoked Motor Response (mV)');
569. %title('Non-Paretic Stimulus-Response Curve'); %label this NP/P as required
570. title(side); %NP/P
571. legend('show','Location','northwest'); % determining the location of the legend
572. box off; % removing the box around the graph
574. %As all tables are stored in one array, here we breaks them back into 4 separate variables
575. FDI = Target_muscles_table{1}; ADM = Target_muscles_table{2}; ECR = Target_muscles_table{3}; FCR = Target_muscles_table{4};
577. % Combining the MEP data from the S-R curve into a table
578. MEPamplitude = table;
579. MEPamplitude.Intensities = Intensities;
580. MEPamplitude.FDI = FDI.MEPamp
581. MEPamplitude.ADM = ADM.MEPamp
582. MEPamplitude.ECR = ECR.MEPamp
583. MEPamplitude.FCR = FCR.MEPamp
584. MEPsum = zeros(1,4);
585. MEPadjusted = MEPamplitude;
586. for muscle = 1:4;
587. subvalue = MEPadjusted(1,muscle+1)
588. m = table2array(MEPadjusted(:,muscle+1))
589. m = m-table2array(subvalue(1,1))
590. MEPadjusted(:,muscle+1) = array2table(m)
591. MEPsum(1,muscle) = sum(table2array(MEPadjusted(:,muscle+1)))
592. end
594. clear adm_bgampstd bgampav colors curveX curveY dy eb ecr_bgampstd fcr_bgampstd fdi_bgampstd fitParams ii initialParams newRowAmp newRowStD sigmoid ADM_temp ECR_temp FCR_temp FDI_temp alpha CI_lower CI_lower_BGamp CI_upper CI_upper_BGamp color_rgb hex_color margin_of_error mean_rms_for_state mean_RMS_state1 MEPamp MEPampArray muscle_name muscle_temp_table muscle_temp_table_name n rows_of_interest s SI std_rms_for_state std_RMS_state1 t_score valid_rows
595. %% Calculating other stats and put them into on struct
596. %this section take the data from the tables and calculate k, meany, df,
597. %SEy, TSS, Critical T and CI. It groups these statistics together with the
598. %corresponding table (for each muscle), storing them in structs.
600. FDI_All = []; ADM_All = []; ECR_All = []; FCR_All = []; %create structs
601. All = {FDI_All, ADM_All, ECR_All, FCR_All}; %store structs in one cell for looping
602. variables = {FDI, ADM, ECR, FCR};
604. for i = 1:length(variables)
605. current_variable = variables{i};
606. current_all = All{i};
607. current_all.Data = current_variable;
608. current_all.SI = SIarray(i,:);
609. current_all.MEPvalues = MEPamplitude(:,i+1);
610. current_all.SigX = SigXarray(i,:); % X values
611. current_all.SigY = SigYarray(i,:); % Y values from the sigmoid fit
612. current_all.S50 = S50(1,i); % assigning S50
613. current_all.k = k(1,i); % k parameter is the 4th value in the parameter output from the sigmoid function
614. current_all.Slope_at_S50 =SlopeAtS50(1,i)*10; %assigning slope at maximum in mV per 10% MSO
616. All{i} = current_all; %assigning the modified struct back into the cell
617. end
619. %breaks the cell back into 4 separate variables
620. FDI = All{1};
621. ADM = All{2};
622. ECR = All{3};
623. FCR = All{4};
625. % Adding the variables you want saved to the Output file
626. Output.FDI = FDI; % All FDI related info
627. Output.ADM = ADM; % All ADM related info
628. Output.ECR = ECR; % All ECR related info
629. Output.FCR = FCR; % All FCR related info
630. Output.MEPamplitude = MEPamplitude; % MEPamp table across intensities and muscles
631. Output.StD = StD; % StDev table across intensities and muscles
632. Output.BGamp = BGamp; % BGamp table across intensities and muscles
633. Output.RMS = RMSmean; % RMS table across intensities and muscles
634. Output.RMSlowerCI = RMSlowerCI; % RMS confidence interval table across intensities and muscles
635. Output.RMSupperCI = RMSupperCI; % RMS confidence interval table across intensities and muscles
636. Output.AverageRMS = avg_rms_values; % average rms per muscle
637. Output.RMS_thresholds = rms_thresholds.'; % the lowesr rms cut off (in 5 microvolt increments) value required for each muscle to include at least 5 frames per intensity
638. Output.k = k.' % adding the k values to the output
639. Output.S50 = S50.'; % S50 values for the four muscles
640. Output.SlopeAtS50_mV1MSO = SlopeAtS50.'; % Slope at S50 in mV/1%MSO for the 4 muscles
641. Output.SlopeAtS50_mV10MSO = SlopeAtS50*10.'; % Slope at S50 in mV/10%MSO for the 4 muscles
642. Output.GoodnessofFit = GoodnessofFit.'; % Goodness of fit of the SR curve for the 4 muscles
643. Output.FrameNumber = frames_per_intensity; % Number of frames per intensity per muscle included for analysis
644. Output.TotalTrials = total_trials;
645. Output.MEPsum = MEPsum;
646. Output.Adjusted_MEPamplitude = MEPadjusted;
647. Output.estimated_threshold = estimated_threshold;
648. Output.adj_estimated_threshold = adj_estimated_threshold;
649. Output.new_adj_estimated_threshold = new_adj_estimated_threshold;
650. Output.percent_retained = percent_retained;
652. clear ADM_All baseline_time current_variable ECR_All FCR_All FDI_All i
653. %% Saving the output file and SR curve figure
655. % Get the current working directory (where the script is located)
656. scriptFolder = pwd;
658. % Define the subfolder for output files
659. outputSubfolder = strcat(assessor, '_Output Files')
661. % Create the full path for the output subfolder
662. outputFolder = fullfile(scriptFolder, outputSubfolder);
664. % Make sure the output subfolder exists
665. if ~exist(outputFolder, 'dir')
666. mkdir(outputFolder);
667. end
669. % Save the analysis Output
670. filename = strcat(fileprefix, 'output', '.mat');
671. fullpath = fullfile(outputFolder, filename);
672. save(fullpath, 'Output');
674. % Save the figure
675. figureFilename = strcat(fileprefix, '_SR Curve.png');
676. figureFullPath = fullfile(outputFolder, figureFilename);
677. saveas(gcf, figureFullPath);
679. % % clearvars -except Output SlopeAtS50 MEPsum

ESPRESSo_SRcurve_current.m at commit e50e60b, no license · at the source