% This model has three species in each finite element
% d_i concentration of the free drug in the ith element
% r_i concentration of the receptor in the ith element
% b_i concentration of the bound receptor in the ith element
speciesTemplates = {'d%d', 'r%d', 'b%d'};
pNames           = {'kon', 'koff', 'D'};
outputTemplates  = {'bound%d'};

% Number of finite elements
noElements       = 50;

% There are three species per element
xNames = cell(noElements*length(speciesTemplates), 1);

% There is one output per element
yNames = cell(noElements,1);

% Build the name strings.
for i = 1:noElements
    for j = 1:length(speciesTemplates)
        xNames{(i-1)*3+j} = sprintf(speciesTemplates{j}, i);
    end
    yNames{i} = sprintf(outputTemplates{1}, i);
end


% Initialize the model structure.  This creates all of the matrices
% of the correct size
m = initModel(length(xNames), 1, noElements,  length(pNames));

m.name   = 'my pde example';

m.xNames = xNames;
m.yNames = yNames;
m.pNames = pNames;

% Output is bound drug in each cell
m.c      = kron(eye(length(xNames)/3),[0,0,1]); 
m.ic(2:3:m.nX) = 1; % Start out with just receptor in each cell

% The system parameters are the on rate, off rate, and diffusion constant.
% The diffusion constant is scaled by the number of elements squared.  This 
% arises because the discretization is made finer, the length of each element
% becomes smaller, and the diffusion constant scales with the square of
% length.
m.p      = [10;1; 0.001*noElements^2];

% Indices of the various variables
drug     = 1:3:m.nX;
receptor = 2:3:m.nX;
complex  = 3:3:m.nX;
kOn      = 1;
kOff     = 2;
D        = 3;

% Connect the model to the input
% drug diffuses into the blood;
m        = addInputOne2OneRev(m, 1, drug(1), D, D);  

% Add the chemistry for each element, and the diffusion between elements
for i = 1:noElements
    
    % In each cell, D + R <--> B
    m = addTwo2OneRev(m, drug(i), receptor(i), complex(i), kOn, kOff);    
    % You can add more reactions here.
    
    if(i ~= noElements)
       % Free drug can diffuse to the next cell
        m = addOne2OneRev(m, drug(i), drug(i+1), D, D);
        % If other things can diffuse, add them here.
    end
end

% postProcessModel finalizes the model, using the reactions defined above.
% This function makes the various matrices that are required for
% model analysis.
m   = postProcessModel(m);



%% Simulate a constant dose of drug.
inputValue = 1;
exp = constructExperiment(m, 1e3, inputValue);
[t,y,x]= simulate(m, exp);



%% Plot bound drug
% In this example the concentration of bound drug is a function of
% time (t) and position (x, or i).  So we will plot a movie of the 
% concentration as a function of position.  Time is indicated on the
% top of the plot.  The subplot on the right shows the same information
% plotted as a color bar
xPlot = linspace(0, 1, m.nY);
for tPlot = linspace(0,t(end), 1000);

    subplot(1,2,1);
    plot(xPlot,y(tPlot));
    axis([0, 1, 0, 1]);
    xlabel('Position (x)');
    ylabel('Bound Receptor [C]');
    title(sprintf('Time = %f\n', tPlot));
    
    
    subplot(1,2,2);    
    image(y(tPlot)*100);
    title(sprintf('Time = %f\n', tPlot));
    drawnow;
end