% NOTES: Created and Copyrighted by William Myles, PhD Candidate at the University of
% Newcastle, Australia. Cite: 'Analytical Methods for Assessing Retinal
% Cell Coupling Using Cut-Loading' Myles and McFadden, 2022.
% If you have any questions or feedback, please don't hesitate to contact
% me, my email is will.myles@newcastle.edu.au 

close all
clear all

%% Data Input
% Mean distance between coupled cells (micrometers)
cellspacing = 62.5;

% Incubation time (seconds)
T = 60;

% Add x coordinates of cells measured (units here are in micrometers)
x = [18.275
34.182
41.3
44.796
102.867
106.315
134.201
160.031
186.07
193.243
235.129
242.542
250.786
254.824
302.684
308.109
356.544
382.76
]; 
% Add corresponding normalised fluorescence data for each cell. Here, the background
% fluorescnece was meaured and subtracted from each data point. Cell fluorescence
% is then presented as a fraction of the maximum fluorescence of cells aligning 
% the cut. These cells aligning the cut are not included in the fit. 

ydata = [0.458012576
0.657280253
0.572215025
0.586812695
0.548906662
0.406187676
0.330962113
0.352412257
0.276622265
0.30018397
0.199710603
0.184267278
0.130699135
0.150582879
0.06046107
0.075700754
0.022540437
0.027234741


];

%% Coupling Coefficient and Diffusion Coefficient Calculation. 
%Depending on the nature of your data, some fit parameters may need to be
%adjusted

%X data normalised to represent number of cells 
xdata = x/cellspacing; 

%X data presented as cm 
xcm = x/10000;

%Calculate De (cm^2.cm^-1)
[Defitresult, Degof] = createDeFit(xcm, ydata);
DeRateConstant = coeffvalues(Defitresult)';

%The equation used for the fit is for tissue incubated for 1500 seconds. As
%the curve fitting app in matlab doesn't let you input values other than
%x,y and z, I have had to adjust the De based on T input after fitting the
%equation. If you like, you can change the time in the associated
%'createDeFit' file. Simply substitute '1500' with your incubation time in
%seconds (you may also have to change the limits of the fitting values). 
%Then simply remove '*1500)/T' from the line below. 
DeRateConstant(4) = (DeRateConstant(2)*1500)/T;
DeRateConstant(3) = DeRateConstant(4)/((cellspacing/10000).^2);

%Create fit to estimate fluorescence at 45 descrete cell units extending from cut
[fitresult gof] = createFitextended(xdata, ydata);
n = 45;
xlin = linspace(1, n, n);
concentrationnew =(fitresult(xlin));

%Calculate kj (cells^2.s^-1) and (cm^2.s^-1)
[Rateconstant, tv] = differentialrateequation(T,cellspacing,concentrationnew);
fitvalues = tv(100,[1:n])';

%Print results for diffusion coefficient fit
fprintf(1,'\tDiffusion Rate Constants:\n')
fprintf(1, '\t\tDe(cells^2.s^-1) = %8.5f\n', DeRateConstant(3));
fprintf(1, '\t\tDe(cm^2.s^-1) = %8.8f\n', DeRateConstant(4));

%Create figure demonstrating fit
figure('Name','Concentrationval');
    plot(xlin,fitvalues);
    set(gca, 'YScale', 'log','fontsize', 14);
    xlim([0 20]);
    ylim([0.0001 1]);
hold on 
    plot(xlin,concentrationnew);
hold on
    scatter(xdata,ydata);
hold off

%%Fit functions 

function [fitresult, gof] = createFitextended(x,y)
[xData, yData] = prepareCurveData(x,y);

% Set up fittype and options.
ft = fittype( 'gauss1' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.StartPoint = [0.02 10 5];
opts.Lower = [0.01 -1 1];
opts.Upper = [2 1 50];

% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );

% Plot fit with data.
figure( 'Name', 'untitled fit 1' );
h = plot( fitresult, xData, yData );
legend( h, 'nonlogyval vs. xdata', 'untitled fit 1', 'Location', 'NorthEast', 'Interpreter', 'none' );
% Label axes
xlabel( 'xdata', 'Interpreter', 'none' );
ylabel( 'nonlogyval', 'Interpreter', 'none' );
grid on
end

function [theta, Cfit] = differentialrateequation(T,CS,IC)
% 22 09 2021
% NOTES: Created and Copyrighted by William Myles, PhD Candidate at the University of
% Newcastle, Australia. Cite: 'Analytical Methods for Assessing Retinal
% Cell Coupling Using Cut-Loading' Myles and McFadden, 2022.


% Arrange data input. Matlab requires three points to fit. An arbitraryily
% small timepoint of 0.0000001s has been included in addition to t = 0 and
% t = T (user defined timepoint when data was collected)
t = [0, 0.0000001, T]';
c = zeros (3,45);
c(:,1) = 8;
c(3,:) = IC;

%Define ode45 solver
function C = kinetics(theta,t)
    c0 = zeros (1,45);
    c0(1,1) = 8;
    [~,Cv]=ode45(@DiffEq,t,c0);

% Define the partial differential equation series. Each PDE models the
% movement of tracer through a subsequent cell in the series. The inclusion
% of theta(3) allows for the initial concentration in C1 to be fixed for
% all fits. This code contains 45 levels, which is typically enough for
% tissues which have been incubated for 20-30 minutes. However, this is
% excessive for tissues with sort incubation times. 
    function dC = DiffEq(~,c)
    dcdt=zeros(5,1);
    dcdt(1)= 2.*theta(1).*(c(2)-c(1))-theta(3).*(c(1));
    dcdt(2)= 2.*theta(1).*(c(1)+c(3)-(2*c(2)))-theta(3).*(c(2));
    dcdt(3)= 2.*theta(1).*(c(2)+c(4)-(2*c(3)))-theta(3).*(c(3));
    dcdt(4)= 2.*theta(1).*(c(3)+c(5)-(2*c(4)))-theta(3).*(c(4));
    dcdt(5)= 2.*theta(1).*(c(4)+c(6)-(2*c(5)))-theta(3).*(c(5));
    dcdt(6)= 2.*theta(1).*(c(5)+c(7)-(2*c(6)))-theta(3).*(c(6));
    dcdt(7)= 2.*theta(1).*(c(6)+c(8)-(2*c(7)))-theta(3).*(c(7)); 
    dcdt(8)= 2.*theta(1).*(c(7)+c(9)-(2*c(8)))-theta(3).*(c(8));
    dcdt(9)= 2.*theta(1).*(c(8)+c(10)-(2*c(9)))-theta(3).*(c(9));
    dcdt(10)= 2.*theta(1).*(c(9)+c(11)-(2*c(10)))-theta(3).*(c(10));
    dcdt(11)= 2.*theta(1).*(c(10)+c(12)-(2*c(11)))-theta(3).*(c(11));
    dcdt(12)= 2.*theta(1).*(c(11)+c(13)-(2*c(12)))-theta(3).*(c(12));
    dcdt(13)= 2.*theta(1).*(c(12)+c(14)-(2*c(13)))-theta(3).*(c(13)); 
    dcdt(14)= 2.*theta(1).*(c(13)+c(15)-(2*c(14)))-theta(3).*(c(14));
    dcdt(15)= 2.*theta(1).*(c(14)+c(16)-(2*c(15)))-theta(3).*(c(15)); 
    dcdt(16)= 2.*theta(1).*(c(15)+c(17)-(2*c(16)))-theta(3).*(c(16));
    dcdt(17)= 2.*theta(1).*(c(16)+c(18)-(2*c(17)))-theta(3).*(c(17)); 
    dcdt(18)= 2.*theta(1).*(c(17)+c(19)-(2*c(18)))-theta(3).*(c(18));
    dcdt(19)= 2.*theta(1).*(c(18)+c(20)-(2*c(19)))-theta(3).*(c(19)); 
    dcdt(20)= 2.*theta(1).*(c(19)+c(21)-(2*c(20)))-theta(3).*(c(20));
    dcdt(21)= 2.*theta(1).*(c(20)+c(22)-(2*c(21)))-theta(3).*(c(21));
    dcdt(22)= 2.*theta(1).*(c(21)+c(23)-(2*c(22)))-theta(3).*(c(22));
    dcdt(23)= 2.*theta(1).*(c(22)+c(24)-(2*c(23)))-theta(3).*(c(23)); 
    dcdt(24)= 2.*theta(1).*(c(23)+c(25)-(2*c(24)))-theta(3).*(c(24));
    dcdt(25)= 2.*theta(1).*(c(24)+c(26)-(2*c(25)))-theta(3).*(c(25)); 
    dcdt(26)= 2.*theta(1).*(c(25)+c(27)-(2*c(26)))-theta(3).*(c(26));
    dcdt(27)= 2.*theta(1).*(c(26)+c(28)-(2*c(27)))-theta(3).*(c(27)); 
    dcdt(28)= 2.*theta(1).*(c(27)+c(29)-(2*c(28)))-theta(3).*(c(28));
    dcdt(29)= 2.*theta(1).*(c(28)+c(30)-(2*c(29)))-theta(3).*(c(29)); 
    dcdt(30)= 2.*theta(1).*(c(29)+c(31)-(2*c(30)))-theta(3).*(c(30));
    dcdt(31)= 2.*theta(1).*(c(30)+c(32)-(2*c(31)))-theta(3).*(c(31));
    dcdt(32)= 2.*theta(1).*(c(31)+c(33)-(2*c(32)))-theta(3).*(c(32));
    dcdt(33)= 2.*theta(1).*(c(32)+c(34)-(2*c(33)))-theta(3).*(c(33)); 
    dcdt(34)= 2.*theta(1).*(c(33)+c(35)-(2*c(34)))-theta(3).*(c(34));
    dcdt(35)= 2.*theta(1).*(c(34)+c(36)-(2*c(35)))-theta(3).*(c(35)); 
    dcdt(36)= 2.*theta(1).*(c(35)+c(37)-(2*c(36)))-theta(3).*(c(36));
    dcdt(37)= 2.*theta(1).*(c(36)+c(38)-(2*c(37)))-theta(3).*(c(37)); 
    dcdt(38)= 2.*theta(1).*(c(37)+c(39)-(2*c(38)))-theta(3).*(c(38));
    dcdt(39)= 2.*theta(1).*(c(38)+c(40)-(2*c(39)))-theta(3).*(c(39)); 
    dcdt(40)= 2.*theta(1).*(c(39)+c(41)-(2*c(40)))-theta(3).*(c(40));
    dcdt(41)= 2.*theta(1).*(c(40)+c(42)-(2*c(41)))-theta(3).*(c(41));
    dcdt(42)= 2.*theta(1).*(c(41)+c(43)-(2*c(42)))-theta(3).*(c(42));
    dcdt(43)= 2.*theta(1).*(c(42)+c(44)-(2*c(43)))-theta(3).*(c(43)); 
    dcdt(44)= 2.*theta(1).*(c(43)+c(45)-(2*c(44)))-theta(3).*(c(44));
    dcdt(45)= 2.*theta(1).*(c(44)+c(45))-theta(3).*(c(45))-theta(2).*(c(45));
    dC=dcdt;
    end
C = Cv;
end

% Set initial estimate for theta 1, 2 and 3
theta0=[0.01;0.001;0.001];

% Fit data with ODE45 solver
[theta,~,~,~,~,~,~]=lsqcurvefit(@kinetics,theta0,t,c);
Kj = (theta(1)*2)*((CS/10000).^2);

% Print rate constant output
fprintf(1,'\tCoupling Rate Constants:\n')
fprintf(1, '\t\tkj(cells^2.s^-1) = %8.5f\n', theta(1));
fprintf(1, '\t\tkj(cm^2.s^-1) = %8.8f\n', Kj(1))


tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);
end

function [fitresult, gof] = createDeFit(xcm, ydata)
%CREATEFIT2(XCM,YDATA)
%  Create a fit.
%
%  Data for 'De_Calculation' fit:
%      X Input : xcm
%      Y Output: ydata
%  Output:
%      fitresult : a fit object representing the fit.
%      gof : structure with goodness-of fit info.
%
%  See also FIT, CFIT, SFIT.

%  Auto-generated by MATLAB on 19-Jan-2022 10:36:02
%% Fit: 'De_Calculation'.
[xData, yData] = prepareCurveData(xcm, ydata);

% Set up fittype and options.
% You can modify this equation to better suit your data. The '0.005' value
% indicates the thickness of l, or the thickness of the cut. For my tissues
% this was typically around 0.005cm or 50um. The '1500' value is the
% incubation time in s. Feel free to substitute your own incubation times
% into the equation. See the notes in the main program for more details. 
ft = fittype( '0.5*C*(erf((0.005-abs(x))/(2*sqrt(D*1500)))+(erf((0.005+abs(x))/(2*sqrt(D*1500)))))', 'independent', 'x', 'dependent', 'y' );
    opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
    opts.Display = 'Off';
    opts.Lower = [0 0];
    opts.StartPoint = [0.913375856139019 0.63235924622541];
    opts.Upper = [100 0.0001];

% Fit De model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );

% Plot fit with data.
figure( 'Name', 'De_Calculation' );
    h = plot( fitresult, xData, yData );
    legend( h, 'ydata vs. xcm', 'De_Calculation', 'Location', 'NorthEast' );
% Label axes
    xlabel xcm
    ylabel ydata
    grid on

end