README.txt for the SPP modeling package as part of the manuscript ÒIdentifying the mechanism for superdiffusivity in mouse fibroblast motilityÓ by Passucci et. al. 

The purpose of this package is to determine if a heterogeneous self-propelled particle model (SPP) or a Levy walk are appropriate models for systems of cells with super diffusive trajectories. The overall procedure relies on using cell trajectory data, so we have provided the ACTIVE nuclei tracking package which identifies and tracks nuclei in stained images to generate trajectories. The ACTIVE package generates the xyzs_id matrix containing the important information of x and y positions, frame numbers and cell identification numbers. 

If you are NOT using ACTIVE, make sure to organize your trajectory data into a four column matrix, where the first column is x positions, the second is y positions, third is frames numbers and the fourth is cell identification numbers. If you are inputting your data in this manner then make sure to set xyzs_id_columns to 4. 

Example Matlab code you could use to generate trajectories: 

dt = 1;
v0 = 1;
a = 1;
ts = 100;
np = 100;
dtmax = ts-1;
posvec = cell(np,1);
for i = 1:np
    rng('shuffle')
    %disp(num2str(i))
    xpos = 0;
    ypos = 0;
    theta = rand(1)*2*pi;
    thetaRan = a * (randn(1,ts));
    posvec{i} = zeros(ts,3);
    for t=1:ts
        dtheta = thetaRan(t);
        theta = theta + dtheta*sqrt(dt);
        xpos = xpos + v0*cos(theta)*dt;
        ypos = ypos + v0*sin(theta)*dt;
        theta = mod(theta,2*pi);
        posvec{i}(t,1) = xpos;
        posvec{i}(t,2) = ypos;
        posvec{i}(t,3) = theta;
    end
end

After appropriately preparing your trajectory data, run data_gather.m with either the output from the ACTIVE code or your own trajectory matrix, along with an estimate of the value of mu, which corresponds to the 3-alpha, where alpha is the slope of the mean-squared displacement on a log-log plot and the first indicator of superdiffusivity. This will generate datas, a matrix containing all of the speed, rotational diffusion, and run times for the run-time segments of your trajectories. You can use this information as a starting location for the parameter sweeps we will talk about later. data_gather.m will also create a file called Òtraj_info.matÓ in your current directory. It is important to be in the same directory as traj_info.mat when you run spp_gen.m. If not then you must modify the code to have an appropriate direct path to that information. This is excluded from the code to be more general since direct paths differ between operating systems. 

Now that you have this basic trajectory information, we can proceed to spp_gen.m. This function creates either a Levy walk (Type = 1) or heterogeneous SPP (Type = 2) with necessary parameters found in vector P, loads traj_info.mat and finds the goodness of fit between the generated statistics and your experimental data. Plot_toggle is set to 0 as a default but if set to any other number spp_gen will produce plots of the aforementioned statistics. I highly suggest experimenting with this function several times before moving on to sweeping phase space. Importantly, the calculations of mean-squared displacement and velocity auto-correlation functions are vectorized, dramatically increasing their efficiency in Matlab. 

A special note on the distributions of parameters in the heterogeneous SPP model. The function mvnRND.m generates a bivariate distribution with correlations of speeds from a 2D Maxwell distribution centered around speed mu_v and normally distributed rotational diffusion constants. If you would like more information on how this was generated please look at the Matlab documentation on creating arbitrary distributions with correlations: https://www.mathworks.com/help/stats/examples/simulating-dependent-random-variables-using-copulas.html#zmw57dd0e233. 

After you are comfortable with the previous steps, proceed to mult_spp_gen.m. Edit the values of p1, p2É. pN appropriately to sweep the desired areas of phase space, where p1É vary depending on which model you are using, determined by the Type parameter. This function will produce figures on how the parameters change around the best goodness-of-fit location in phase space as well as save the important information in a file in the same directory that this function is run. 

Finally, we have prt_collapse, which collapses displacement probability distributions using pcell and collapse_times (generated by data_gather.m for experimental data and spp_gen.m for simulation data). An important input parameter for this function is ntrials, as this controls the number of iterations for a bootstrap estimate of error propagation on the probability distributions. Distributions are constructed using the ksdensity Matlab function which we found efficient and accurate. 

Acknowledgements: I am in no way responsible for the code contained in the ACTIVE tracking package, otsu.m, or the edge_detector_1d package. The vectorization of MSD calculations was inspired by the function Kohl.m. At this time of this writing these external functions are available either the Matlab file exchange or at the request of the research groups responsible for them.