Reaction-Diffusion Model and Data VisualizationJustin MitchellBlockedUnblockFollowFollowingFeb 24In the scientific community, a model is a representation of a concept that can be used to simulate and help describe an idea.
If the model reflects the real world closely enough, it can also be used to make predictions about what would happen under a set of conditions.
One representation of a model that is commonly encountered is a weather map, as seen below.
National Centers for Environmental Prediction [Public domain], via Wikimedia CommonsA weather map is the visual representation of a weather model.
A weather model takes inputs such as temperature, pressure, wind speed, and other parameters at many locations and, after processing this information, it can be used to describe and predict the weather.
Another aspect of models to consider is that they can be wrong.
Anyone who has foregone taking an umbrella as they leave the house based on a weather prediction knows this.
Therefore, it is important to keep in mind that models are often imperfect recreations of the real world.
However, this does not mean that they can not be useful.
For example, when a physical system is well understood, a model can be developed to allow experiments within the system to be run in a virtual environment.
This allows scientists to explore the variation of parameters of the system, which may allow them to discover an optimal set of conditions for a given target.
The reason one might opt to create a model rather than physically run experiments on a system is that it can oftentimes be expensive and time consuming to gather enough data over a range of parameters to be cost effective in the task of exploring optimal system conditions.
One such system is a reaction diffusion model.
A chemical system that involves reaction and diffusion has many variables and when the components involved are costly, it can be infeasible to run trials over the entire parameter domain.
A model of this type of system will be presented in this post.
Specifically, the model in this post will represent the reaction and diffusion occurring in a cubic chemical well with active sites adsorbed to the walls and base of the well.
The wells are similar to the ones pictured below, but the geometry will be cubic instead of cylindrical.
The surfaces of the well are filled with reaction sites that the chemicals in the well solution interact with.
Obtained from Wikimedia CommonsObtained from https://www.
comThe reaction space, a cube, will be divided into smaller cubes in the same manner as the cube pictured to the left.
The dimensions of the large cube will be 10x10x10 small cubes, hereafter referred to as cells.
The model will work by developing an array with four dimensions: x, y, z, and t.
Respectively, these represent the x, y, and z spatial dimensions, and the time dimension.
The values at each array index will represent the concentration at that space and time.
The two important factors in a reaction diffusion model are, unsurprisingly, reaction and diffusion.
The reaction takes place only at the walls of the well but diffusion occurs throughout the system.
After reacting, the chemical is considered to be removed from the system.
After diffusing, the chemical moves from one well to another, but is not removed from the system.
First the system must be defined with initial conditions.
These include a diffusivity constant, the length of each cell side, length of simulation time, length of time step, initial concentration of the cells, initial well surface site concentration, the rate constants of adsorption and desorption, and initial concentration of adsorbed chemical.
At time 0, all cells are equal to the initial concentration.
The cell concentrations for the next time step are calculated by removing the chemical lost due to reaction from each cell and calculating the changes in cell concentration from diffusion based on the difference in concentrations between the cell of interest and all neighboring cells.
The python code completing this task is presented below.
"""D = Diffusivity constant in m^2/sr = length of each dimension of cubic well in metersn = number of partitions in each dimension of cubic welltf = length of time of simulationdt = numerical integration time stepinitial_conc = initial concentraion of well solutiondx = r/nnt = tf/dtxrn is in moles per timeAconc and Sconc are in moles per m^2k_des is in per timek_ad is in per M per tdx is in mrxnDMtot is total change in cube molarity due to reaction"""stracker = 0D = 1 * 10 ** -5r = 0.
0048tf = 0.
1dt = 0.
000001n = 10nt = int(tf / dt)dx = r/ninitial_conc = 2.
5 * 10 ** -5import numpy as npconc = np.
zeros([n,n,n,nt+5])for i in range(0, n): for j in range(0, n): for k in range(0, n): conc[i, j, k, 0] += initial_conc#print(conc)Sconc = 1 * 10 ** -8 #surface site concentrationAconc = np.
zeros([nt+6]) #adsorbed concentrationDfactor = D *dt/(dx ** 2)k_ad = 1 * 10 ** 7k_des = 1 * 10 ** -6Sconc = 1Aconc = 0time_list = Aconc_list = SminusA_conc_list = for t in range(0, nt + 6): for i in range(0, n): for j in range(0,n): for k in range(0,n): #start of neighbor unit diffusion calculations if i == n – 1: xposdiff = 0 else: xposdiff = Dfactor * (conc[i + 1, j, k, t] – conc[i, j, k, t])if i == 0: xnegdiff = 0 else: xnegdiff = Dfactor * (conc[i – 1, j, k, t] – conc[i, j, k, t]) if j == n – 1: yposdiff = 0 else: yposdiff = Dfactor * (conc[i, j + 1, k, t] – conc[i, k, k, t]) if j == 0: ynegdiff = 0 else: ynegdiff = Dfactor * (conc[i, j – 1, k, t] – conc[i, j, k, t]) if k == n – 1: zposdiff = 0 else: zposdiff = Dfactor * (conc[i, j, k + 1, t] – conc[i, j, k, t]) if k == 0: znegdiff = 0 else: znegdiff = Dfactor * (conc[i, j, k – 1, t] – conc[i, j, k, t]) #end of neighbor unit diffusion calculations #start of neighbor unit reaction calculations if i == 0 or i == n – 1: xrxn = dx ** 2 * (k_ad * conc[i, j, k, t] * (Sconc – Aconc[t]) – k_des * Aconc[t]) #xrn is in moles per time, Aconc and Sconc are in moles per m^2, k_des is in per time, k_ad is in per M per t, dx is in m else: xrxn = 0 if j == 0 or j == n – 1: yrxn = dx ** 2 * (k_ad * conc[i, j, k, t] * (Sconc – Aconc[t]) – k_des * Aconc[t]) else: yrxn = 0if k == 0: zrxn = dx ** 2 * (k_ad * conc[i, j, k, t] * (Sconc – Aconc[t]) – k_des * Aconc[t]) else: zrxn = 0 #end of neighbor unit reaction calculations#calculates total molarity change due to reactions rxnDMtot = -(xrxn + yrxn + zrxn) * dt / (dx * 10) ** 3 #total change in cube molarity due to reaction#limits change in concentration due to reaction because time step minimization is limited by processing power if rxnDMtot < -1 * conc[i, j, k, t]: rxnDMtot = -1 * conc[i, j, k, t] #keeps track of the surface site concentration stracker = stracker – rxnDMtot * (dx * 10) ** 3#uses the current unit concentration and the reaction and diffusion changes to calculate the unit concentration at the next time step conc[i, j, k, t + 1] = conc[i, j, k, t] + xposdiff + xnegdiff + yposdiff + ynegdiff + zposdiff + znegdiff + rxnDMtot Aconc[t + 1] = Aconc[t] + stracker / (5 * r ** 2) stracker = 0 time_list.
append(t * dt) Aconc_list.
append(Sconc – Aconc[t])for i in range(0, n): for j in range(0,n): for k in range(0,n): totalconc = totalconc + conc[i, j, k, nt] * (dx * 10) ** 3totalconc = totalconc / ((r * 10) ** 3)totalconc = totalconc + Aconc[nt] * 5 * r ** 2 / (r * 10) ** 3One of the main reasons I undertook this project was because I wanted to develop my skills in presenting information using python.
Displaying the four dimensional array that contains cell concentrations would require five dimensions to display the three spatial dimensions, concentrations, and time.
This could be accomplished with an animated three dimensional heat map.
I may return to this project to accomplish displaying all the information at a later time.
For now, I have developed gifs using equally-spaced-in-time snapshots of the cell concentrations with locked z dimensions.
What this means is that the x and y axis will represent the location of the cells in the x and y direction, the vertical axis will represent the concentration of each cell, and the time dimension will be represented by the animation.
The z coordinate being displayed is indicated under each gif.
The second from the bottom layer (z = 1)The fourth from the base (z = 3)You can see that the concentration in the cells at the sides of the reaction space drops to zero immediately.
We’ll refer to these cells as the wall cells since they are in contact with the wall of the reaction space.
The reason the concentrations in the wall cells drop faster than the interior cells is because the wall cells are the only cells that are losing concentration due to reaction.
Diffusion results in solute flow from high concentrations to low concentrations.
Due to being higher in concentration than the wall cells, diffusion causes the concentrations of the interior cells to begin to drop after reaction begins because solute has a net transfer toward the walls.
Unfortunately, this is not a perfect model because there are much more sophisticated methods to describe this system more accurately and the time step is not small enough to accurately reflect reality.
AcknowledgementsJohn Roberts for his work in developing the original VBA model.
Luke Borsare for his assistance in the development of the graphics.
Fang Fang, Kevin, and Caroline for their guidance.
Graphics and Gif Generation Python Codeimport matplotlib!pip install imageioimport mpl_toolkitsfrom mpl_toolkits.
mplot3d import axes3dimport matplotlib.
pyplot as plt!pip install –upgrade matplotlibdef plot_looper(snapshot_number, z_level): #this will only work correctly with n = 10 t = 0 file_title = 'z' while t < tf / dt: fig = plt.
figure() ax = fig.
add_subplot(111, projection = '3d') X = [[0,1,2,3,4,5,6,7,8,9],[9,0,1,2,3,4,5,6,7,8],[8,9,0,1,2,3,4,5,6,7],[7,8,9,0,1,2,3,4,5,6],[6,7,8,9,0,1,2,3,4,5],[5,6,7,8,9,0,1,2,3,4],[4,5,6,7,8,9,0,1,2,3],[3,4,5,6,7,8,9,0,1,2],[2,3,4,5,6,7,8,9,0,1],[1,2,3,4,5,6,7,8,9,0]] Y = [[0,1,2,3,4,5,6,7,8,9],[0,1,2,3,4,5,6,7,8,9],[0,1,2,3,4,5,6,7,8,9],[0,1,2,3,4,5,6,7,8,9],[0,1,2,3,4,5,6,7,8,9],[0,1,2,3,4,5,6,7,8,9],[0,1,2,3,4,5,6,7,8,9],[0,1,2,3,4,5,6,7,8,9],[0,1,2,3,4,5,6,7,8,9],[0,1,2,3,4,5,6,7,8,9]] Z = conc[X,Y,z_level, t] file_title += 'z' ax.
scatter(X, Y, Z, c = 'r', marker = 'o') fig.
format('PATH/x')) print(t) t += 1000 import osimport imageiopng_dir = 'PATH/Model_Images'images = for file_name in os.
listdir(png_dir): if file_name.
png'): file_path = os.
join(png_dir, file_name) images.
gif', images).. More details