lipowskiModel/analyseSimDynamics.m at main · HoldenLab/lipowskiModel · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
function [simResult, chgpt, tFr,xFr] = analyseSimDynamics(t,x,frameInterval,zeroSpeedThresh, switchTooCloseThresh, runID,dynamics)
% function [simResult, chgpt, tFr,xFr] = analyseSimDynamics(t,x,frameInterval,zeroSpeedThresh, switchTooCloseThresh)
%INPUTS
% frame interval: units: t (typically s): exposure time for a single image frame in comparable experimental system. For MreB dynamics, we usually use 6 s
% zeroSpeedThresh - units: x/t (typically nm/s): speed below which speed cannot be accurately differentiated from stationary molecule
% switchTooCloseThresh - units: frames, number of frames below which a switching event is too close to be detected in the comparable experimental analysis. Typically 4 frames based on things like MSD analysis

if ~exist('runID','var')
    runID=1;
end

frameInterval=frameInterval;
vThresh=zeroSpeedThresh;
frameThresh = switchTooCloseThresh;

tFr = 0:frameInterval:max(t);
xFr = interp1(t,x,tFr);

dxFr= diff(xFr);
v=dxFr/frameInterval;

%filter out near zero speed
v2= v;
v2(abs(v2)<=vThresh)=0;

%loop through and find sign changes or pauses
%need to identify states(beginning and end)
%measure type of state (motile, paused), length of state, average speed, and next state
%(motile paused)
%exclude states that are too close to last one
vSign = v2;
vSign(vSign>0) =1;
vSign(vSign<0) =-1;

%plot(tFr(1:end-1),vSign);
chgpt=find(diff(vSign)~=0);
%add the first and last frames as states
chgpt=[1,chgpt,numel(tFr)-1];

tChgPt=tFr(chgpt);
%plot(tChgPt,vSign(chgpt),'x')
%Filter out short trajectories
chgptFilt=chgpt;
ii = 2;
while ii <= numel(chgptFilt)
    if ii > 1 && (chgptFilt(ii)-chgptFilt(ii-1))<frameThresh
        if chgptFilt(ii)+1>numel(vSign) %if last change point out of bounds
            %delete last chgptFilt
            chgptFilt(ii)=[];
        elseif vSign(chgptFilt(ii-1))==vSign(chgptFilt(ii)+1)%previous state same as current state
            %delete both chgptFilt
            chgptFilt(ii-1:ii)=[];
            ii=ii-1;
        else
            %dlete only previous chgptFilt
            chgptFilt(ii-1)=[];
        end
    else
        ii=ii+1;
    end
end

%add back in the first and last points if they are sufficiently
%far away from a change point
if chgptFilt(1)~=1 && chgptFilt(1)-1>=frameThresh
    chgptFilt=[1,chgptFilt];
end
if chgptFilt(end)~=numel(vSign) && (numel(vSign)-chgptFilt(end))>=frameThresh
    chgptFilt=[chgptFilt, numel(vSign)];
end

if exist('dynamics','var')
    NplusFr = interp1(t,dynamics.Nplus,tFr);
    NminusFr = interp1(t,dynamics.Nminus,tFr);
    NmotorFr = NplusFr+NminusFr;
end


%record states in output table
simResult = table;
nState = numel(chgptFilt)-1;
simResult.runID = ones(nState,1)*runID;
%loop through and measure state properties
for ii=1:nState
    startIdx = chgptFilt(ii);
    endIdx  = chgptFilt(ii+1);
    simResult.duration(ii)= tFr(endIdx)-tFr(startIdx);
    simResult.processivity(ii)=abs(xFr(endIdx)-xFr(startIdx));
    simResult.speed(ii)=simResult.processivity(ii)/simResult.duration(ii);
    simResult.isMotile(ii)=vSign(endIdx)~=0;
    if exist('dynamics','var')
        simResult.Nmotor(ii) = median(NmotorFr(startIdx:endIdx));
    end
    if ii<nState
        nextStateIdx=chgptFilt(ii+2);
        simResult.isNextStateMotile(ii)=double(vSign(nextStateIdx)~=0);
    else
        simResult.isNextStateMotile(ii)=NaN;
    end
end

chgpt = chgptFilt;