
% This program creates an XML file in the current directory. The file can
% be run in OpenCell to simulate the neuron model given on p.400 in Butera      
% et. al 1999b. g_syn_e is given pseudo-random values from a normal
% distribution of user specified mean and standard deviation.
%
% Note: To run the file in OpenCell, the one cell model from the repository
% must be in the same folder, and its name should not have been changed.

clc ; clear ;

% Ask user how many neurons are to be modelled
disp (' ') ; neuronCount = input ('How many neurons are to be modelled : ') ;

% Prompt user to enter the mean and standard deviation of the normal dist.
disp (' ') ;
g_syn_e_mean = input ('Please enter the mean of g_syn_e normal distribution : ') ;
g_syn_e_stdDev = input ('Please enter the standard deviation of g_syn_e normal distribution : ') ;

% Create array of pseudo-random g_syn_e values
g_syn_e = g_syn_e_mean + g_syn_e_stdDev * randn (neuronCount, neuronCount) ;

% Begin the XML code
xmlCode = sprintf ('<?xml version="1.0"?> \n<model xmlns="http://www.cellml.org/cellml/1.1#" xmlns:cmeta="http://www.cellml.org/metadata/1.0#" cmeta:id="butera_%d_cell_1999" name="butera_%d_cell_1999"> \n\n', neuronCount, neuronCount) ;

% Adding imports and connections 
for i = 1 : neuronCount
    imports = sprintf ('  <import xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="butera_single_cell_1999.cellml"> \n    <component component_ref="single_neuron_model" name="single_neuron_model%d"/> \n  </import> \n\n  <connection> \n    <map_components component_1="single_neuron_model%d" component_2="synaptic_coupling"/> \n    <map_variables variable_1="s" variable_2="s%d"/> \n    <map_variables variable_1="sum_g_syn_e_s" variable_2="sum_g_syn_e_s_%d"/> \n  </connection> \n\n', i, i, i,i) ;
    xmlCode = [xmlCode imports] ;
end

% Adding g_syn_e units and initialising the synaptic_coupling component code
g_syn_e_stuff = sprintf ('  <units name="nanoS"> \n    <unit prefix="nano" units="siemens"/> \n  </units> \n\n  <component cmeta:id="synaptic_coupling" name="synaptic_coupling"> \n') ;
xmlCode = [xmlCode g_syn_e_stuff] ;

% Adding sum_g_syn_e_s variable descriptions
for i = 1 : neuronCount
    sum_g_syn_e_s_descriptors = sprintf ('    <variable name="sum_g_syn_e_s_%d" public_interface="out" units="nanoS"/> \n', i) ;
    xmlCode = [xmlCode sum_g_syn_e_s_descriptors] ;
end

% Formatting
space = sprintf ('\n') ;
xmlCode = [xmlCode space] ;

% Adding g_syn_e values
for i = 1 : neuronCount
    for j = 1 : neuronCount
        if (i ~= j)
            g_syn_e_values = sprintf ('    <variable initial_value="%f" name="g_syn_e_%d_%d" units="nanoS"/> \n', g_syn_e (i, j), i, j) ;
            xmlCode = [xmlCode g_syn_e_values] ;
        end
    end
end

% Formatting
space = sprintf ('\n') ;
xmlCode = [xmlCode space] ;

% Adding g_syn_e components
for i = 1 : neuronCount
    g_syn_e_components = sprintf ('    <variable name="s%d" public_interface="in" units="dimensionless"/> \n', i) ;
    xmlCode = [xmlCode g_syn_e_components] ;
end

% Formatting
space = sprintf ('\n') ;
xmlCode = [xmlCode space] ;

% Initialising the code for g_syn_e math
initialMath = sprintf ('    <math xmlns="http://www.w3.org/1998/Math/MathML"> \n') ;
xmlCode = [xmlCode initialMath] ;

% Adding g_syn_e math code
for j = 1 : neuronCount
    sum_g_syn_e_s_math = sprintf ('      <apply><eq /> \n        <ci> sum_g_syn_e_s_%d </ci> \n        <apply><plus /> \n', j) ;
    xmlCode = [xmlCode sum_g_syn_e_s_math] ;
    for i = 1 : neuronCount
        if (i ~= j)
            g_syn_e_math = sprintf ('		  <apply><times/> \n            <ci> g_syn_e_%d_%d </ci> \n            <ci> s%d </ci> \n          </apply> \n', i, j, i) ;
            xmlCode = [xmlCode g_syn_e_math] ;
        end
    end
    mathEnd = sprintf ('		</apply> \n	  </apply> \n\n') ;
    xmlCode = [xmlCode mathEnd] ;
end

% Ending the g_syn_e component code
g_syn_e_end = sprintf ('	</math> \n  </component> \n\n</model>') ;
xmlCode = [xmlCode g_syn_e_end] ;

% Create XML file
myfid = fopen ('buteraModelXML.xml','w') ;

% Write code to file
fprintf (myfid, '%s', xmlCode) ;

% Close file
fclose(myfid) ;
