Simple Cell Traction Force Script for Elastic Micropatterned Substrata

Hi readers, 

 This is MATLAB code that was written by myself and collaborators that we've sought to make available to the wider research community. The program is intended to track the displacements of micropatterned dots on a substrate in a similar manner to that performed by  Maloney et al. in "Influence of Finite Thickness on cellular adhesion-induced deformation of an compliant substrata". Physical Review E. 2008.

 The program is only intended to be for MATLAB and uses programs from other members of the research/MATLAB community to track the displacements of labeled "spots" on a surface of a substrate, sufficiently spaced (approx 5um or greater) so that they do not interfere with eachother. We hope to eventually incorporate a number of other details, such as the application of forces that will influence neighboring points or for points that are more closely spaced than 5um and to increase the level of automation the program has.

Thanks for reading and trying it out! Please leave suggestions if you have any. We will post a link to the article, demonstrating the technique and methods on how to apply it once it's published. 

 

%%Script for Traction Force Calculations Made on an Elastic Gel%%%

%%Version 5.0 - Written and Contributions by: Samuel Polio, Katheryn Rothenberg, Dimitrije Stamenovic, and Michael Smith

%%Boston University, Deparment of Biomedical Engineering

%%

 

%uncomment for clearing everything before calculating 

% close all; 

% clear; 

m = [];

pk = []; 

cnt = [];

a = [];

b = [];

pos = [];

 

%Note that values need to be entered with YOUR parameters, please adjust accordingly where the XX's are 

 mcircle_dia = XX; %maximum circle diameter in um

 post_dist = XX; %distance between posts in um

 pixel_ratio = XX; %number of um per pixel

pixel_dia = round(mcircle_dia / pixel_ratio); %max diameter of pixels for circle

pix_post_dist = post_dist/pixel_ratio; %distance between posts in pixels

 

%Get the dot file!

initial_file = uigetfile('*.tif','Pick The Dot File'); 

 

%read in .tif file

initial_image = double(imread(initial_file, 'tif')); 

%display the image

colormap(gray);

imagesc(initial_image);

 

%please freezeColors from: 

%http://www.mathworks.com/matlabcentral/fileexchange/7943-freezecolors-unfreezecolors

freezeColors

 

%georgetown macros... please obtain from http://physics.georgetown.edu/matlab/

    bp = bpass(initial_image,5,pixel_dia); %bandpass image to remove background

    m = max(max(bp))*.1; % Percent of the brightest feature estimate

    pk = pkfnd(bp, m, pixel_dia); %find all peak locations

    cnt = cntrd(bp,pk,pixel_dia+2); %accurately locate peak centroid

 

hold;

plot(cnt(:,1),cnt(:,2),'x');

 

figure 

plot(cnt(:,1),cnt(:,2),'x'); 

 

 

%selecting points on plot for the corners

datacursormode on; 

fig_handle = gcf;

hold on; 

pos = zeros(4,2);

for i=1:4,

    choice = menu('Press Enter After Selecting Point','Enter');

    while choice==0

        choice = menu('Press Enter After Selecting Point','Enter');

    end

%     pause(); 

    dcm_obj = datacursormode(fig_handle); 

    info = getCursorInfo(dcm_obj);

    pos(i,1) = info.Position(1);

    pos(i,2) = info.Position(2); 

    plot(pos(i,1),pos(i,2),'rx'); 

end

 

%putting the 4 corners into the program... Find this program below 

[x,y]=FourCorners(pos(1,:),pos(2,:),pos(3,:),pos(4,:));

 

%setting time to 0 and 1 so that the program knows which dots to connect

xy = [x,y]; 

xy(:,3) = 1;

 

cnt(:,3)=[]; 

cnt(:,3)=0; 

 

t_pos = [cnt;xy]; 

tracked = track(t_pos, 30); 

 

count = 0; 

for i = 1: (length(tracked)-1)

    if tracked(i,4) == tracked (i+1,4)

        count = count+1; 

        p1(count,1) = tracked(i,1);

        p1(count,2) = tracked(i,2);

        p2(count,1) = tracked(i+1,1);

        p2(count,2) = tracked(i+1,2);

    end

end

 

%"Error" is if there is no displacement, but one can concieve of it to be the

% magnitudes of the displacements

error = sqrt((p1(:,1)-p2(:,1)).^2+(p1(:,2)-p2(:,2)).^2);

 

cell_image_file = uigetfile('*.tif','Pick The Dot File', 'F:\03_30_11'); %plotting image for the cell with arrows

 

%display the figure

figure; 

colormap('gray')

%uncomment ONE of the following at a time to show the forces

%for displaying the dot displacement, uncomment 

imagesc(double(imread(initial_file)))

%for displaying the displacement on the brightfield image, uncomment

%imagesc(double(imread(initial_file))) 

hold; 

freezeColors

 

A = p1 - p2; %displacements for the image

 

thresh2 = 0; %threshold displacement of dots in pixels

under = find(error < thresh2); %locating the points under the threshold to set to 0

 

for i = 1:length(under), 

    A(under(i),1) = 0; 

    A(under(i),2) = 0; 

end

 

colormap('jet')

%forces are shown at the final position and scaled, if you put the dot at 

%the initial position (p2), then unscale with quiverc(x,y,u1,u2,0), you can show 

%the actual displacements of the dots

%quiverc from Bertrand Dano, http://www.mathworks.com/matlabcentral/fileexchange/3225-quiverc 

quiverc(p1(:,1), p1(:,2), A(:,1), A(:,2)); 

 

%Get force/displacement vectors magnitudes

Disp = error*pixel_ratio; 

 

%Force magnitudes from gel mechanical data

%Calculate the spring constant from Maloney et al. "Influence of Finite Thickness

%on cellular adhesion-induced deformation of compliant substrata". Physical Review E. 2008. 

%For larger dots or dots that interfere with each other even though they 

%are not being pulled on, different parameters need to be used 

nu = XX;%Poisson's Ratio

dot_radius = XX; %meters

Elastic_mod = XX; %pascals

k = pi()*(Elastic_mod*dot_radius)/((1+nu)*(2-nu));%spring constant of gel in N/m

k_nano = k*10^3; %spring constant in nN/um

Force = Disp*k_nano;%nN

 

%Vectors for the force and displacement

D_Vec = A*pixel_ratio; 

F_Vec = D_Vec*k_nano; 

 

%%SETTING UP PLOT AXES%%

N = 5; 

max_F = max(Force); %max total Force and cell number, can change to displacement to show that as well. 

sep = max_F/N; %separation between ticks - (denominator) - 1

 

%Hard-coded axis labels

num1 = num2str(0,'%6.2f'); 

num2 = num2str(sep,'%6.2f'); 

num3 = num2str(sep*2,'%6.2f'); 

num4 = num2str(sep*3,'%6.2f');

num5 = num2str(sep*4,'%6.2f');

num6 = num2str(sep*5,'%6.2f');

 

hcb = colorbar('YTickLabel', {num1,num2,num3,num4,num5,num6});

 

cax_vals = caxis; 

 

cax_width = (cax_vals(2)-cax_vals(1))/5;

 

set(hcb,'YTick', [cax_vals(1),cax_vals(1)+cax_width,cax_vals(1)+2*cax_width...

    cax_vals(1)+3*cax_width,cax_vals(1)+4*cax_width,cax_vals(1)+5*cax_width]); 

%%END PROGRAM %%

 

 

 

 

 

FourCorners.m

 

function [Xf Yf] = FourCorners(pos1, pos2, pos3, pos4, distance)

 

%The function FourCorners takes in 4 sets of xy pairs representing the 4

%corners of a grid of dots. The idea is to use these 4 corners to calibrate

%the distance between the dots in 2 orthogonal axes. Program needs a slight tile to work properly. 

%distance = distance in pixels between dots

 

%Find the corners

P = [pos1; pos2; pos3; pos4]; 

 

dist = zeros(length(P)-1,1); 

i=1;

%finding the distance between each point, the longest distance between 2

%points will be on a diagonal and used to not knock out the pairing between

%the first point and that one when drawing the grid

for i = 2:(length(P))

    j = i-1; 

    dist(j,1) = i;

    dist(j,2) = sqrt((P(1,1)-P(i,1))^2+(P(1,2)-P(i,2))^2);

end

 

[~,ma] = max(dist(:,2)); %on diagonal

[~,mi] = min(dist(:,2)); %value don't want to use for next one

 

%find the longest distance not on the diagonal to get the best average 

%along the length

for i = 1:length(dist), 

    if dist(i,2) ~= max(dist(:,2)) && dist(i,2) ~= min(dist(:,2));  

        conn = dist(i,1); 

    end

end

 

set1 = [pos1;P(conn,:)]; %arrange data into sets for theoretically parallel lines

set2 = [P(dist(mi,1),:); P(dist(ma,1),:)];

 

%get equations for both lines

s1 = polyfit(set1(:,1),set1(:,2),1); 

s2 = polyfit(set2(:,1),set2(:,2),1); 

 

%manual input

%ask user number of dots between the points 

s = sprintf('How many spots are between points %d,%d and %d,%d?\n',...

    set1(1,1), set1(1,2),set1(2,1), set1(2,2));

num_parallel = input(s); 

 

%automatic input

%num_s = sqrt((set1(1,1)-set1(2,1)).^2 +(set1(1,2)-set1(2,2)).^2);

%num_p = (num_s)/distance; 

%num_parallel = round(num_p); 

 

%calculate the x separation between the points on the parallel lines

total_d_parallel_1x = set1(1,1)-set1(2,1);

sep_parallel_1x = abs(total_d_parallel_1x / (num_parallel+1)); 

 

total_d_parallel_2x = set2(1,1)-set2(2,1);

sep_parallel_2x = abs(total_d_parallel_2x / (num_parallel+1)); 

 

 

 

%Preallocating for speed... 

x1 = zeros(num_parallel+2,1); 

y1 = zeros(num_parallel+2,1); 

x2 = zeros(num_parallel+2,1); 

y2 = zeros(num_parallel+2,1); 

 

%"draw" grid points along each parallel line to dot first line

for i = 1:num_parallel+2, 

    if set1(1,1)<set1(2,1) %start a point from the left and go to the right

        x1(i) = set1(1,1)+(i-1)*sep_parallel_1x;

        y1(i) = s1(1)*x1(i)+s1(2); 

    elseif set1(1,1)>set1(1,2)

        %go backwards if it's not right

        x1(i) = set1(1,1)-(i-1)*sep_parallel_1x;

        y1(i) = s1(1)*x1(i)+s1(2);

    end

end

 

%"draw" grid points along each parallel line to dot second line

for i = 1:num_parallel+2, 

    if set2(1,1)<set2(2,1) %start a point from the left and go to the right

        x2(i) = set2(1,1)+(i-1)*sep_parallel_2x;

        y2(i) = s2(1)*x2(i)+s2(2); 

    elseif set2(1,1)>set2(1,2)

        %go backwards if it's not right

        x2(i) = set2(1,1)-(i-1)*sep_parallel_2x;

        y2(i) = s2(1)*x2(i)+s2(2);

    end

end

 

%manual calculate the number of dots to put between

s = sprintf('How many spots are between points %d,%d and %d,%d?\n',...

    set1(1,1), set1(1,2),set2(1,1), set2(1,2));

num_perp = input(s); 

 

%automatic input

%num_perp = round((sqrt((set1(1,1)-set2(1,1))^2 +(set1(1,2)-set2(1,2))^2))/distance); 

 

Xf = []; 

Yf = []; 

 

for i = 1:length(x1)

    s = polyfit([x1(i);x2(i)],[y1(i);y2(i)],1);

    sep_perp_x = abs((x1(i)-x2(i))/(num_perp+1)); 

 

    for k = 1:num_perp+2,

        if x1(i)<x2(i) %start a point from the left and go to the right

            X(k) = x1(i)+(k-1)*sep_perp_x;

            Y(k) = s(1)*X(k)+s(2); 

        elseif x1(i)>x2(i)

        %go backwards if it's not right

            X(k) = x1(i)-(k-1)*sep_perp_x;

            Y(k) = s(1)*X(k)+s(2); 

        end

    end

    

    Xf = [Xf;X']; 

    Yf = [Yf;Y']; 

end

 

end