This is the repository of code associated with the paper, "Concentration-Dependent Domain Evolution in Reaction-Diffusion Systems." It is divided into an implementation of the 1D model in MATLAB, and an implementation of the 1D and 2D model in COMSOL.
The MATLAB code can be used in two ways. Firstly, you can run the main.m script to generate an output.mat structure which will be automatically saved in a local file named output/modelParamsSelector/growthParamsSelector/output.m (here using the Linux file separator, but the directory will be created correctly independently of OS). Alternatively, you can run the sweep.m file after uncommenting out a selection of kinetics and growth parameters to run multiple simultaneous simulations using the MATLAB Parallel Computing Toolbox.
Each simulation is described by two strings called modelParamsSelector and growthParamsSelector. The string modelParamsSelector corresponds to diffusion rates, kinetic functions and parameters, steady states, and initial conditions. The growthParamsSelector corresponds to domain length and growth rate parameters. Starting at main.m, you can look at the corresponding functions where these are specified to see what the defaults are, or to add/modify these just by adding a new 'case' under the switch statements in getModelParams.m and getGrowthParams.m modelled after the existing cases. If the runUniform flag (defined in the growthParams structure) is set to 1, main.m or sweep.m will also generate a simulation using the same beginning and end domain lengths to run over the same timescale but at a uniform exponential rate - that is, with S = const. This will be stored in the file output/modelParamsSelector/growthParamsSelector/uniform/output.m.
The function plotRun(FILEPATH) will create a kymograph of the concentration of the first variable ('u' in the paper) and store it in the same directory as output.m, as well as three snapshots in time similar to those shown in Figure 3 in the paper. The functions writeVideo and animateRun can also be used to create some very interesting videos or live animations of the solutions in order to better understand how they behave.
Some important implementation details: the function S must return a vector, so in all cases where we specify a constant value of S, we also add +0u to achieve this. The simulationParams structure contains mesh size and timestep controls, which can be modified via setupSimulationOptions.m. Initial conditions are by default set as random perturbations of the provided steady state values, but these can be modified by changing the modelParams.initType string in getModelParams.m and then writing a corresponding initial condition in setupModel.m.
These have been provided for the 2D simulations, as well as for comparisons done between the 1D and the 2D codes. Start by looking at one of the directories given by Concentration-Dependent-Growth-Simulations/COMSOL 1D and 2D Code/2D/Simulations/NAME, where NAME is one of the folders corresponding to the examples provided. The ceateDomain.m file produces the outermost boundary curve of the initial simulation domain. From here, using MATLAB LiveLink, we ran the file runModel.m which should automatically step through solving the model through a number of iterations as described in that file, and saving the resulting values of the Eulerian coordinates and u in the NAME/csvs/ folder. These csvs can then be visualized in a number of different ways as described in the files in the plottingCodes folder.