"
Team:TU-Munich/Modeling/Filter
From 2013.igem.org
Filter calculator
To simplify the testing of use cases and for the comfort of future users, we implemented a filter calculator. The filter calculator computes the number of remediation rafts needed to degrade a certain amount of pollutant in different environmental conditions. To do this we characterised the kinetics of the laccase at different pH values, temperatures and in different salt concentrations (see our kinetic fitting section for details). Alternatively you can enter the parameters of your own protein. For more information on how this is computed see the documentation below.
Calculate number of remediation rafts
Select enzyme:
Select water way:
River flow rate [m3/s]:
Medium flow speed [m/s]:
Enter pH of water way:
Enter water temperature [°C]:
Enter salt concentration [M]:
Enter enzyme activity (kcat/KM) [M-1s-1]:
Enter enzyme half-life [h]:
Enter enzyme mass [Da]:
Subtrate concentration before filter [nM]:
Distance of downstream measurement after filter [km]:
Desired substrate concentration after filter in [nM]:
Filter model
Aim
We wanted to know, how our PhyscoFilter could be used in practice and what kind of scale would be needed in a real-world implementation. So we looked at the some example water ways and waste water treatment plants, taking into account not only the local contamination, but also the characteristics of the water, such as pH, salt concentration and temperature, as well as flow speed and volume. To do this, we programmed an interactive calculator computing the required number of rafts for the given setup and provided example data for different use-cases with our laccase.
Model
To model this we assumed that the enzymes are secreted into the flowing water of the river or outflow of a waste water treatment plant. To determine the enzyme concentration we could achieve in this way, we questioned Prof. Reski and the experts of Greenovation (see our Expert council), a company using P. patens for the mass production of proteins for the pharmaceutical market. They told us, that in their experience fully grown moss produces 10-50 mg of protein per liter of liquid culture per day. For our computations we chose the conservative value of 25 mg per liter of culture per day.
Our remediation rafts (see our implementation page) have an area of slightly over 1 m2, so taking into account the additional density due to the sedentary growth on the felt, we estimated this to roughly equal a liquid culture of 10 cm height, giving a culture volume of 100 liters per raft.
The price per raft is (also on our implementation page) about 50 $ including the setup. The optional Arduino controller providing live data on the current status of the filter, costs about 250 $, but you only need a single one for each filter location.
The enzymes produced in a day will be diluted in the amount of water flowing by in a day, so the initial enzyme concentration in the water is given by:
Now it must be taken into account, that the enzymes are inactivated over time, so given an inactivation rate λ or a half-life T1/2, we have:
For the evolution of the substrate concentration we started with the well-known Michaelis-Menten formula and dropped the substrate concentration in the denominator, since the concentrations of the micro-pollutants found in the environment is always much smaller than the KM values of the proteins, so we have the differential equation
Integration gives 
By inverting this we get the formula below, showing how many rafts are necessary to reduce the substrate concentration by a specified amount
Since it is easier to work with the downstream distance along the river than with the time since the water passed the filter, we related the two according to 
.
pH and salt dependence
The inactivation rates of our enzyme in different pH-values and salt concentrations (see our enzyme kinetic fitting) were clearly not realistic, since the inactivation was too fast, and in strong disagreement with our assay of the laccase inactivation in river water. However the trends seemed sensible, so we decided to normalize the inactivation rates to fit with the half-life observed in the river water, which had pH 7.5 and a salt concentration of 0.28 mM (from this water analysis).
Assuming the inactivation rates to be additive, we obtained the plot in figure 2 showing the dependence of the inactivation rate on the pH-value and salt concentration.
River data values
As an example substrate to illustrate the calculator we chose the painkiller Diclofenac, which is widely used and very stable in the environment, hence reaching significant concentrations in the environment. It is tragically famous for being responsible for the near extinction of Indian vultures, in whose bodies high concentrations had accumulated after eating livestock excessively treated with Diclofenac. Diclofenac can be oxidised and hence inactivated by the laccase, so it is a great target for our PhyscoFilter.
The data on the Diclofenac pollution and the general characteristics of the rivers was obtained from: The Rhine in Germany: [BLAC, 2003], [Monitory station Rhine], [RIWA, 2008].
The Kaveri in India: [Shanmugam et al., 2013], [Begum and Harikrishnarai, 2008].
A large WWTP in Germany (Großlappen near Munich): [BLAC, 2003], [Waste water treatment plant Großlappen report], [Report on drinking water in the area Isar-Vils].
Conclusion
As an example for a large waste water treatment plant in Germany the calculator tells us, that about 3000 remediation rafts are needed to reduce the Diclofenac (a painkiller, responsible for the near extinction of vultures in India) from the average value found there to levels below the suggested EU limit. This is quite a lot, but the area required would be less than half a football (i.e. soccer) pitch. For larger water ways, where the volume of the passing water dilutes the enzymes a lot, the calculator returns even higher numbers, which was to be expected. However, the activity of the laccase could be optimized and the literature values on the kinetics of the laccase from Bacillus pumilus vary widely both above and below our values. For example, [Thöny-Meyer et al., 2011] obtained a laccase activity for ABTS of 3.6 * 106, so 20 times our estimate, which would greatly reduce the number of remediation rafts required. The model also illustrates very nicely the different dependencies, which play together in the complex degradation process combined with the enzyme inactivation. A good example the change in the required number of rafts, when the temperature is increased, indicating good applicability in warmer countries, such as India, but also potential problems in countries with high temperature amplitudes. All in all, the model shows, that this approach could work on water ways with higher concentrations and relatively low flow rates, but that some optimization is required to deal with large rivers, where the enzyme is extremely diluted.
Code:
The code for this calculator can be found here.
Diffusion model
Aim
With this model we wanted to get an idea of how the cell wall and the cell membrane would influence the diffusion of a pollutant. Additionally the aim was to investigate the effect of the enzyme localization on the filter efficiency.
Model
We simulated an idealized two-dimensional section of a cell with polluted water flowing past on the outside of the cell and pollutants diffusing inward through the cell wall and cell membrane. Then we added enzymes at the different locations our constructs support, namely in the cytoplasm, bound to the cell membrane or secreted into the intercellular space (assuming the protein is so large it does not pass through the cell wall).
Results
Figure 3 shows the animated evolution of the pollutant concentration in the absence of degrading enzymes and assuming the pollutant freely diffuses through the cell wall, while continuously pollutant is washed in from the left. As expected pollutant accumulates over time and will eventually reach the constant concentration of the in-flow.
Figure 4 and 5 show the simulation for cytoplasmic and intercellular enzyme localization respectively, again assuming free diffusion through the cell wall and keeping all other parameters (including the enzyme concentration) as before. The difference compared to Figure 3 is obvious and it is clearly visible (and intuitive), that the intercellular enzymes have a greater effect.
To show the further options of our program, the simulation of the effect of membrane bound enzymes with the cell wall impeding the diffusion of the pollutant is shown in Figure 6. However, since the enzymes are now spread out over a surface instead of a volume this is not properly comparable to the previous figures.
Conclusion
In general the model is not quantitative at the moment, because parameters like the diffusivity of the cell membrane etc. are very hard to determine. However, the simulations created give a good indication and a nice visualization of what is happening.
The Mathematics behind it
The evolution of the concentration of the pollutant is governed by the Advection-Diffusion-Equation describing the combination of advective, i.e. through movement of the medium, and diffusive transport of a substance in a medium such as water.
where D is the diffusion coefficient, v is the flow speed and R is the degradation rate of the enzymes, proportional to the concentration and the activity of the enzyme.
We approximated the right-hand side of this equation by breaking up the area into a finite grid and using finite differences (as shown here). On the boundaries inside the cell we assumed that there is no concentration flow across, i.e. 
, where the vector n is normal to the boundary. To simulate the exchange across the membrane and where appropriate the cell wall, we used
where k is a constant proportional to the diffusivity of the membrane.
On the left boundary outside of the cell wall, the concentration was kept fixed simulating a constant inflow of pollutant. On the remaining boundaries the diffusive term across the boundary was dropped, since the advective flux is much greater than the diffusive one there.
The resulting differential equation in time was solved using a MatLab ODE solver.
Code
Main script: Diffusion model
filename = 'NoWallNoEnzyme.gif'; % where to save the output
xLength = 80; % must be an integer, but should be greater 20 for proper resolution,
% but program takes approx. 1 minute at xlength=100 and time goes up rapidly
widthOfExtracellularSpace = 20; % must be an integer, same size constraints as above
widthOfIntercellularSpace = 10; % must be an integer
widthOfCytoplasmaticSpace = 20; % must be an integer
flowSpeed = 2; % speed with which fluid flows on outside of cell
diffusionWater = 0.005; % diffusion parameter in water
CellWallEffect = 0; % 0, if the cell wall does NOT impede diffusion at all
% 1, if the cell wall impedes diffusion
diffusionCellWall = 30; % diffusion parameter for the cell wall
% value NOT used, if CellWallEffect = 0, but needs some value for easier programming
diffusionCellMembrane = 10; % diffusion parameter for the cell membrane
TimeHorizon = 5; % total time to be modeled
TimeSteps = 0.2; % step size for the clip
% (TimeHorizon/TimeSteps + 1 gives the total number of frames)
enzymeType = 0; %enzymeType takes the following values: 0 - for no enzymes
% 1 - for cytoplasmic enzymes
% 2 - for enzymes fixed to the membrane
% 3 - for enzymes in the intercellular space
enzymeRate = 100; % degradation rate of the enzymes, including concentration etc.
%=== CALCULATIONS FROM HERE ONWARDS, DO NOT CHANGE ANYTHING BELOW THIS LINE ===
totalWidth = widthOfExtracellularSpace + widthOfIntercellularSpace + widthOfCytoplasmaticSpace;
params = createParamVec(xLength,widthOfExtracellularSpace,widthOfIntercellularSpace,flowSpeed,diffusionWater,diffusionCellWall,diffusionCellMembrane,enzymeType,enzymeRate);
initialSetup = [[ones(widthOfExtracellularSpace,1);zeros(widthOfIntercellularSpace + widthOfCytoplasmaticSpace,1)],zeros(totalWidth,xLength-1)];
disp('Setup complete, beginning calculations. May take a few minutes')
if CellWallEffect == 1
[TOut,YOut] = ode113(@(t,y) timeDeriv(y,params), 0:TimeSteps:TimeHorizon, reshape(initialSetup,totalWidth*xLength,1));
elseif CellWallEffect == 0
[TOut,YOut] = ode113(@(t,y) timeDeriv_FreeWall(y,params), 0:TimeSteps:TimeHorizon, reshape(initialSetup,totalWidth*xLength,1));
else
error('CellWallEffect must be either 0 (for perfect diffusion across wall) or 1 (for diffusion impeded by wall)')
end
disp('Calculations completed, generating output.')
input('Press Enter to produce the output! Do NOT click away till process is done')
for i = 1 : length(TOut)
if (abs(min(YOut(i,:))) < 10^(-5)) % so there might be output in the negatives, but this is very small, so set it to 0 for the plot
YOut(i,(YOut(i,:)<0)) = 0;
end
[~,h] = contourf(flipud(reshape(YOut(i,:),totalWidth,xLength)),0:0.01:1);
set(h,'LineStyle','none'); % to make the output nicer
rectangle('Position',[0,widthOfCytoplasmaticSpace-0.5,xLength,1],'FaceColor','k')
rectangle('Position',[0,totalWidth-widthOfExtracellularSpace-0.5,xLength,1],'FaceColor','k')
text(xLength/3,totalWidth - widthOfExtracellularSpace/2,'Extracellular Domain','FontSize',18);
text(xLength/3,widthOfCytoplasmaticSpace + widthOfIntercellularSpace/2,'Intercellular Domain','FontSize',18);
text(xLength/3,widthOfCytoplasmaticSpace/2,'Cytoplasmic Domain','FontSize',18);
caxis([0 1]) % fix the colour scale
if (enzymeType == 1) % cytoplasmatic degradation enzymes
% draw enzymes in cyto
rectangle('Position',[ xLength/4-2, widthOfCytoplasmaticSpace/3-2,4,4],'Curvature',[1,1],'FaceColor','k');
rectangle('Position',[2*xLength/4-2,2*widthOfCytoplasmaticSpace/3-2,4,4],'Curvature',[1,1],'FaceColor','k');
rectangle('Position',[3*xLength/4-2, widthOfCytoplasmaticSpace/3-2,4,4],'Curvature',[1,1],'FaceColor','k');
rectangle('Position',[ xLength/4-2,2*widthOfCytoplasmaticSpace/3-2,4,4],'Curvature',[1,1],'FaceColor','k');
rectangle('Position',[2*xLength/4-2, widthOfCytoplasmaticSpace/3-2,4,4],'Curvature',[1,1],'FaceColor','k');
rectangle('Position',[3*xLength/4-2,2*widthOfCytoplasmaticSpace/3-2,4,4],'Curvature',[1,1],'FaceColor','k');
end
if enzymeType == 2 % degradation enzymes localized on the cell mmbrane
rectangle('Position',[ xLength/4-2,widthOfCytoplasmaticSpace-1,4,4],'Curvature',[1,1],'FaceColor','k');
rectangle('Position',[2*xLength/4-2,widthOfCytoplasmaticSpace-1,4,4],'Curvature',[1,1],'FaceColor','k');
rectangle('Position',[3*xLength/4-2,widthOfCytoplasmaticSpace-1,4,4],'Curvature',[1,1],'FaceColor','k');
end
if enzymeType == 3 % degradation enzymes transported into the intercellular space
rectangle('Position',[ xLength/4-2,widthOfCytoplasmaticSpace + widthOfIntercellularSpace/3-2,4,4],'Curvature',[1,1],'FaceColor','k');
rectangle('Position',[2*xLength/4-2,widthOfCytoplasmaticSpace + 2*widthOfIntercellularSpace/3-2,4,4],'Curvature',[1,1],'FaceColor','k');
rectangle('Position',[3*xLength/4-2,widthOfCytoplasmaticSpace + widthOfIntercellularSpace/3-2,4,4],'Curvature',[1,1],'FaceColor','k');
rectangle('Position',[ xLength/4-2,widthOfCytoplasmaticSpace + 2*widthOfIntercellularSpace/3-2,4,4],'Curvature',[1,1],'FaceColor','k');
rectangle('Position',[2*xLength/4-2,widthOfCytoplasmaticSpace + widthOfIntercellularSpace/3-2,4,4],'Curvature',[1,1],'FaceColor','k');
rectangle('Position',[3*xLength/4-2,widthOfCytoplasmaticSpace + 2*widthOfIntercellularSpace/3-2,4,4],'Curvature',[1,1],'FaceColor','k');
end
set(gcf,'Color','w');
drawnow
frame = getframe(gcf);
im = frame2im(frame);
[imind,cm] = rgb2ind(im,256);
if i == 1
imwrite(imind,cm,filename,'gif', 'Loopcount',inf);
else
imwrite(imind,cm,filename,'gif','WriteMode','append');
end
end
% movie2avi(F,filename,'fps',5);
Sub function: Time Derivative with no effect of cell wall on pollutant
function [ concChange ] = timeDeriv_FreeWall( concVector, parameterVector )
%timeDeriv_FreeWall - gives the time derivative of the concentration, when
%the wall does NOT impede diffusion
% parameterVector = [100,40,20,1,1,0.5,0.002];
%% Setup
length = parameterVector(1);
N = length; %grid size
totalWidth = size(concVector,1)/length;
outerWidth = parameterVector(2);
interWidth = parameterVector(3);
cytoWidth = totalWidth - outerWidth - interWidth;
u = parameterVector(4); %flow speed
D = parameterVector(5); %diffusion parameter in water
k_wall = parameterVector(6); %diffusion parameter through cell wall
k_membrane = parameterVector(7); %diffusion parameter through cell membrane
%the diffusion parameters across the wall & membrane are given by:
% permeability*surfaceArea/((1m^3)*thickness of the membrane)
%NOTE: surfaceArea DEPENDS ON N, but depends on NORMALIZATION (i.e. length scale of x-axis)
%the extent in the not shown z-direction is taken to be 1m (1 meter)
%then the formula used for the cross membrane/wall diffusion is:
% dc_inside/dt = - dc_outside/dt = ( permeability*surfaceArea/((1m^3) * thickness) ) * (c_outside - c_inside)
% = k_membrane * (c_outside - c_inside)
enzymeType = parameterVector(8); %0 for none, 1 for cytoplasmic, 2 for localized, 3 for intercellular space
enzymeRate = parameterVector(9); %reaction rate for the enzymes (including concentration)
%transform to matrix for better handling
c = reshape(concVector, totalWidth , length); %c the matrix of concentrations
cOuter = c( 1 : outerWidth , :);
cInter = c( outerWidth + 1 : outerWidth + interWidth , :);
cCyto = c(totalWidth - cytoWidth +1 : totalWidth , :);
% cNew = zeros(totalWidth, length); %the new matrix for storing the output
%% Cores of the domains (i.e. non-boundary terms)
%OUTER DOMAIN
newOuterCore = D*(N^2)*(...
cOuter(3:end,2:end-1) + cOuter(1:end-2,2:end-1) + cOuter(2:end-1,3:end) + cOuter(2:end-1,1:end-2) ...
-4 * cOuter(2:end-1,2:end-1)... %approximation of Laplacian(c)
)...
- (u*N/2)*(...
cOuter(2:end-1, 3:end) - cOuter(2:end-1, 1:end-2)... %approximation of dc/dx
);
%INTER DOMAIN
newInterCore = D*(N^2)*(...
cInter(3:end,2:end-1) + cInter(1:end-2,2:end-1) + cInter(2:end-1,3:end) + cInter(2:end-1,1:end-2) ...
-4 * cInter(2:end-1,2:end-1)... %approximation of Laplacian(c)
);
%CYTO DOMAIN
newCytoCore = D*(N^2)*(...
cCyto(3:end,2:end-1) + cCyto(1:end-2,2:end-1) + cCyto(2:end-1,3:end) + cCyto(2:end-1,1:end-2) ...
-4 * cCyto(2:end-1,2:end-1)... %approximation of Laplacian(c)
);
%% Boundary terms
%OUTER DOMAIN
newOuterLeft = zeros(outerWidth,1); %since it is constant = 1
newOuterRightCore = D*(N^2)*(...
-2*cOuter(2:end-1,end) + cOuter(1:end-2,end) + cOuter(3:end,end))...
-(u*N/2)*(...
3*cOuter(2:end-1,end) - 4*cOuter(2:end-1, end-1) + cOuter(2:end-1,end-2) ...
);
newOuterTopRight = -(u*N/2)*(3*cOuter(1,end) - 4*cOuter(1,end-1) + cOuter(1,end-2));
newOuterTopCore = D*(N^2)*(...
- 2*cOuter(1,2:end-1) + cOuter(1,1:end-2) + cOuter(1,3:end)...
)...
- (u*N/2)*(...
cOuter(1,3:end) - cOuter(1,1:end-2)...
);
newOuterBottomRight = D*(N^2)* (...
cOuter(end-1,end) -2*cOuter(end,end) + cInter(1,end)...
)...
-(u*N/2)*(3*cOuter(end,end) - 4*cOuter(end,end-1) + cOuter(end,end-2)) ...
;
newOuterBottomCore = D*(N^2)*(...
cOuter(end-1,2:end-1) - 4*cOuter(end,2:end-1) + cOuter(end,1:end-2) + cOuter(end,3:end) + cInter(1,2:end-1)...
)...
- (u*N/2) * ( cOuter(end,3:end) - cOuter(end,1:end-2) ) ...
;
%INTER DOMAIN
newInterLeftCore = D*(N^2)*(...
-4*cInter(2:end-1,1) + 2*cInter(2:end-1,2) + cInter(1:end-2,1) + cInter(3:end,1) );
newInterRightCore = D*(N^2)*(...
-4*cInter(2:end-1,end) + 2*cInter(2:end-1,end-1) + cInter(1:end-2,end) + cInter(3:end,end) );
newInterTopLeft = D*(N^2)* (...
-4*cInter(1,1) + 2*cInter(1,2) + cInter(2,1) + cOuter(end,1) ...
);
newInterTopRight = D*(N^2)* (...
-4*cInter(1,end) + 2*cInter(1,end-1) + cInter(2,end) + cOuter(end,end)...
);
newInterTopCore = D*(N^2)*(...
-4*cInter(1,2:end-1) + cInter(2,2:end-1) + cInter(1,1:end-2) + cInter(1,3:end) + cOuter(end,2:end-1)...
);
newInterBottomRight= D*(N^2)* (...
-4*cInter(end,end) + 2*cInter(end,end-1) + 2*cInter(end-1,end) ...
) + ( k_membrane * (cCyto(1,end) - cInter(end,end)) ) ... % diffusion across the membrane
;
newInterBottomCore = D*(N^2)*(...
-4*cInter(end,2:end-1) + 2*cInter(end-1,2:end-1) + cInter(end,1:end-2) + cInter(end,3:end) ...
) + ( k_membrane * (cCyto(1,2:end-1) - cInter(end,2:end-1)) )... % diffusion across the membrane
;
newInterBottomLeft = D*(N^2)* (...
-4*cInter(end,1) + 2*cInter(end,2) + 2*cInter(end-1,1) ...
) + ( k_membrane * (cCyto(1,end) - cInter(end,end)) ) ... % diffusion across the membrane
;
%CYTO DOMAIN
newCytoLeftCore = D*(N^2)*(...
-4*cCyto(2:end-1,1) + 2*cCyto(2:end-1,2) + cCyto(1:end-2,1) + cCyto(3:end,1) );
newCytoRightCore = D*(N^2)*(...
-4*cCyto(2:end-1,end) + 2*cCyto(2:end-1,end-1) + cCyto(1:end-2,end) + cCyto(3:end,end) ...
);
newCytoTopLeft = D*(N^2)* (...
-4*cCyto(1,1) + 2*cCyto(1,2) + 2*cCyto(2,1) ...
) + ( k_membrane * (cInter(end,1) - cCyto(1,1)) ) ... % diffusion across the membrane
;
newCytoTopRight = D*(N^2)* (...
-4*cCyto(1,end) + 2*cCyto(1,end-1) + 2*cCyto(2,end)...
) + ( k_membrane * (cInter(end,end) - cCyto(1,end)) ) ... % diffusion across the membrane
;
newCytoTopCore = D*(N^2)*(...
-4*cCyto(1,2:end-1) + 2*cCyto(2,2:end-1) + cCyto(1,1:end-2) + cCyto(1,3:end) ...
) + ( k_membrane * (cInter(end,2:end-1) - cCyto(1,2:end-1)) ) ... % diffusion across the membrane
;
newCytoBottomRight= D*(N^2)*(...
-4*cCyto(end,end) + 2*cCyto(end-1,end) + 2*cCyto(end,end-1) ...
);
newCytoBottomCore = D*(N^2)*(...
-4*cCyto(end,2:end-1) + 2*cCyto(end-1,2:end-1) + cCyto(end,1:end-2) + cCyto(end,3:end) ...
);
newCytoBottomLeft = D*(N^2)*(...
-4*cCyto(end,1) + 2*cCyto(end-1,1) + 2*cCyto(end,2) ...
);
%% Put together cNew
cNewOuter = [ newOuterLeft , [newOuterTopCore;newOuterCore;newOuterBottomCore] , [newOuterTopRight;newOuterRightCore;newOuterBottomRight] ];
cNewInter = [[newInterTopLeft;newInterLeftCore;newInterBottomLeft] , [newInterTopCore;newInterCore;newInterBottomCore] , [newInterTopRight;newInterRightCore;newInterBottomRight] ];
cNewCyto = [[ newCytoTopLeft; newCytoLeftCore; newCytoBottomLeft] , [ newCytoTopCore; newCytoCore; newCytoBottomCore] , [ newCytoTopRight; newCytoRightCore; newCytoBottomRight] ];
%% add enzymatic activity
if enzymeType == 0
cNewFull = [cNewOuter ; cNewInter ; cNewCyto];
elseif enzymeType == 1 % cytoplasmatic degradation enzymes
cNewFull = [cNewOuter ; cNewInter ; cNewCyto - enzymeRate*cCyto];
elseif enzymeType == 2 % degradation enzymes localized on the cell membrane
cNewInter(end,:) = cNewInter(end,:) - enzymeRate*cInter(end,:);
cNewFull = [cNewOuter ; cNewInter ; cNewCyto];
elseif enzymeType == 3 % degradation enzymes transported into the intercellular domain
cNewFull = [cNewOuter ; cNewInter - enzymeRate*cInter ; cNewCyto];
else
error('enzymeType must be one of: 0 - for no enzymes /n 1 - for cytoplasmic enzymes /n 2 - for enzymes fixed to the membrane /n 3 - for enzymes in the intercellular space')
end
%% transform back to vector for output
concChange = reshape(cNewFull, size(concVector,1),1);
end
Sub function: Time Derivative with the cell wall as diffusion barrier for the pollutant
function [ concChange ] = timeDeriv( concVector, parameterVector )
%timeDeriv computes the time derivative of the concentration, when the cell
%wall impedes diffusion
% concVector must be a COLUMN vector
%% Setup
length = parameterVector(1);
N = length; %grid size
totalWidth = size(concVector,1)/length;
outerWidth = parameterVector(2);
interWidth = parameterVector(3);
cytoWidth = totalWidth - outerWidth - interWidth;
u = parameterVector(4); %flow speed
D = parameterVector(5); %diffusion parameter in water
k_wall = parameterVector(6); %diffusion parameter through cell wall
k_membrane = parameterVector(7); %diffusion parameter through cell membrane
%the diffusion parameters across the wall & membrane are given by:
% permeability*surfaceArea/((1m^3)*thickness of the membrane)
%NOTE: surfaceArea DEPENDS ON N, but depends on NORMALIZATION (i.e. length scale of x-axis)
%the extent in the not shown z-direction is taken to be 1m (1 meter)
%then the formula used for the cross membrane/wall diffusion is:
% dc_inside/dt = - dc_outside/dt = ( permeability*surfaceArea/((1m^3) * thickness) ) * (c_outside - c_inside)
% = k_membrane * (c_outside - c_inside)
enzymeType = parameterVector(8); %0 for none, 1 for cytoplasmic, 2 for localized, 3 for intercellular space
enzymeRate = parameterVector(9); %reaction rate for the enzymes (including concentration)
%transform to matrix for better handling
c = reshape(concVector, totalWidth , length); %c the matrix of concentrations
cOuter = c( 1 : outerWidth , :);
cInter = c( outerWidth + 1 : outerWidth + interWidth , :);
cCyto = c(totalWidth - cytoWidth +1 : totalWidth , :);
% cNew = zeros(totalWidth, length); %the new matrix for storing the output
%% Cores of the domains (i.e. non-boundary terms)
%OUTER DOMAIN
newOuterCore = D*(N^2)*(...
cOuter(3:end,2:end-1) + cOuter(1:end-2,2:end-1) + cOuter(2:end-1,3:end) + cOuter(2:end-1,1:end-2) ...
-4 * cOuter(2:end-1,2:end-1)... %approximation of Laplacian(c)
)...
- (u*N/2)*(...
cOuter(2:end-1, 3:end) - cOuter(2:end-1, 1:end-2)... %approximation of dc/dx
);
%INTER DOMAIN
newInterCore = D*(N^2)*(...
cInter(3:end,2:end-1) + cInter(1:end-2,2:end-1) + cInter(2:end-1,3:end) + cInter(2:end-1,1:end-2) ...
-4 * cInter(2:end-1,2:end-1)... %approximation of Laplacian(c)
);
%CYTO DOMAIN
newCytoCore = D*(N^2)*(...
cCyto(3:end,2:end-1) + cCyto(1:end-2,2:end-1) + cCyto(2:end-1,3:end) + cCyto(2:end-1,1:end-2) ...
-4 * cCyto(2:end-1,2:end-1)... %approximation of Laplacian(c)
);
%% Boundary terms
%OUTER DOMAIN
newOuterLeft = zeros(outerWidth,1); %since it is constant = 1
newOuterRightCore = D*(N^2)*(...
-2*cOuter(2:end-1,end) + cOuter(1:end-2,end) + cOuter(3:end,end))...
-(u*N/2)*(...
3*cOuter(2:end-1,end) - 4*cOuter(2:end-1, end-1) + cOuter(2:end-1,end-2) ...
);
newOuterTopRight = -(u*N/2)*(3*cOuter(1,end) - 4*cOuter(1,end-1) + cOuter(1,end-2));
newOuterTopCore = D*(N^2)*(...
- 2*cOuter(1,2:end-1) + cOuter(1,1:end-2) + cOuter(1,3:end)...
)...
- (u*N/2)*(...
cOuter(1,3:end) - cOuter(1,1:end-2)...
);
newOuterBottomRight = D*(N^2)* (... %here use dc/dy = 0 across membrane, since diffusion across is covered below
-2*cOuter(end,end) + 2*cOuter(end-1,end)...
)...
-(u*N/2)*(3*cOuter(end,end) - 4*cOuter(end,end-1) + cOuter(end,end-2)) ...
+ ( k_wall * (cInter(1,end) - cOuter(end,end)) ) ... % diffusion across the wall
;
newOuterBottomCore = D*(N^2)*(...
2*cOuter(end-1,2:end-1) - 4*cOuter(end,2:end-1) + cOuter(end,1:end-2) + cOuter(end,3:end) ...
) - (u*N/2) * ( cOuter(end,3:end) - cOuter(end,1:end-2) ) ...
+ ( k_wall * (cInter(1,2:end-1) - cOuter(end,2:end-1)) ) ... % diffusion across the wall
;
%INTER DOMAIN
newInterLeftCore = D*(N^2)*(...
-4*cInter(2:end-1,1) + 2*cInter(2:end-1,2) + cInter(1:end-2,1) + cInter(3:end,1) );
newInterRightCore = D*(N^2)*(...
-4*cInter(2:end-1,end) + 2*cInter(2:end-1,end-1) + cInter(1:end-2,end) + cInter(3:end,end) );
newInterTopLeft = D*(N^2)* (...
-4*cInter(1,1) + 2*cInter(1,2) + 2*cInter(2,1) ...
) + (k_wall * (cOuter(end,1) - cInter(1,1)) ) ... % diffusion across the wall
;
newInterTopRight = D*(N^2)* (...
-4*cInter(1,end) + 2*cInter(1,end-1) + 2*cInter(2,end) ...
) + ( k_wall * (cOuter(end,end) - cInter(1,end)) ) ... % diffusion across the wall
;
newInterTopCore = D*(N^2)*(...
-4*cInter(1,2:end-1) + 2*cInter(2,2:end-1) + cInter(1,1:end-2) + cInter(1,3:end) ...
) + ( k_wall * (cOuter(end,2:end-1) - cInter(1,2:end-1)) ) ... % diffusion across the wall
;
newInterBottomRight= D*(N^2)* (...
-4*cInter(end,end) + 2*cInter(end,end-1) + 2*cInter(end-1,end) ...
) + ( k_membrane * (cCyto(1,end) - cInter(end,end)) ) ... % diffusion across the membrane
;
newInterBottomCore = D*(N^2)*(...
-4*cInter(end,2:end-1) + 2*cInter(end-1,2:end-1) + cInter(end,1:end-2) + cInter(end,3:end) ...
) + ( k_membrane * (cCyto(1,2:end-1) - cInter(end,2:end-1)) )... % diffusion across the membrane
;
newInterBottomLeft = D*(N^2)* (...
-4*cInter(end,1) + 2*cInter(end,2) + 2*cInter(end-1,1) ...
) + ( k_membrane * (cCyto(1,end) - cInter(end,end)) ) ... % diffusion across the membrane
;
%CYTO DOMAIN
newCytoLeftCore = D*(N^2)*(...
-4*cCyto(2:end-1,1) + 2*cCyto(2:end-1,2) + cCyto(1:end-2,1) + cCyto(3:end,1) );
newCytoRightCore = D*(N^2)*(...
-4*cCyto(2:end-1,end) + 2*cCyto(2:end-1,end-1) + cCyto(1:end-2,end) + cCyto(3:end,end) ...
);
newCytoTopLeft = D*(N^2)* (...
-4*cCyto(1,1) + 2*cCyto(1,2) + 2*cCyto(2,1) ...
) + ( k_membrane * (cInter(end,1) - cCyto(1,1)) ) ... % diffusion across the membrane
;
newCytoTopRight = D*(N^2)* (...
-4*cCyto(1,end) + 2*cCyto(1,end-1) + 2*cCyto(2,end)...
) + ( k_membrane * (cInter(end,end) - cCyto(1,end)) ) ... % diffusion across the membrane
;
newCytoTopCore = D*(N^2)*(...
-4*cCyto(1,2:end-1) + 2*cCyto(2,2:end-1) + cCyto(1,1:end-2) + cCyto(1,3:end) ...
) + ( k_membrane * (cInter(end,2:end-1) - cCyto(1,2:end-1)) ) ... % diffusion across the membrane
;
newCytoBottomRight= D*(N^2)*(...
-4*cCyto(end,end) + 2*cCyto(end-1,end) + 2*cCyto(end,end-1) ...
);
newCytoBottomCore = D*(N^2)*(...
-4*cCyto(end,2:end-1) + 2*cCyto(end-1,2:end-1) + cCyto(end,1:end-2) + cCyto(end,3:end) ...
);
newCytoBottomLeft = D*(N^2)*(...
-4*cCyto(end,1) + 2*cCyto(end-1,1) + 2*cCyto(end,2) ...
);
%% Put together cNew
cNewOuter = [ newOuterLeft , [newOuterTopCore;newOuterCore;newOuterBottomCore] , [newOuterTopRight;newOuterRightCore;newOuterBottomRight] ];
cNewInter = [[newInterTopLeft;newInterLeftCore;newInterBottomLeft] , [newInterTopCore;newInterCore;newInterBottomCore] , [newInterTopRight;newInterRightCore;newInterBottomRight] ];
cNewCyto = [[ newCytoTopLeft; newCytoLeftCore; newCytoBottomLeft] , [ newCytoTopCore; newCytoCore; newCytoBottomCore] , [ newCytoTopRight; newCytoRightCore; newCytoBottomRight] ];
%% add enzymatic activity
if enzymeType == 0
cNewFull = [cNewOuter ; cNewInter ; cNewCyto];
elseif enzymeType == 1 % cytoplasmatic degradation enzymes
cNewFull = [cNewOuter ; cNewInter ; cNewCyto - enzymeRate*cCyto];
elseif enzymeType == 2 % degradation enzymes localized on the cell mmbrane
cNewInter(end,:) = cNewInter(end,:) - enzymeRate*cInter(end,:);
cNewFull = [cNewOuter ; cNewInter ; cNewCyto];
elseif enzymeType == 3 % degradation enzymes transported into the intercellular space
cNewFull = [cNewOuter ; cNewInter - enzymeRate*cInter ; cNewCyto];
else
error('enzymeType must be one of: 0 - for no enzymes /n 1 - for cytoplasmic enzymes /n 2 - for enzymes fixed to the membrane /n 3 - for enzymes in the intercellular space')
end
%% transform back to vector for output
concChange = reshape(cNewFull, size(concVector,1),1);
end
References:
[BLAC, 2003] Report on pharmaceuticals in the environment
[Monitory station Rhine] Monitory station Karlsruhe for the river Rhine
[RIWA, 2008] Report on pollution in the Rhine basin by RIWA Netherlands
[Shanmugam et al., 2013] Non-steroidal anti-inflammatory drugs in Indian rivers. Environmental Science and Pollution Research, July 2013
[Abida Begum and Harikrishnarai, 2008] Begum, A. and Harikrishnarai (2008). Study on the Quality of Water in Some Streams of Cauvery River. E-Journal of Chemistry, 5 (2008), Issue 2: 377-384.
[Waste water treatment plant Großlappen report] Report about the Waste water treatment plant Großlappen
[Report on drinking water in the area Isar-Vils]
[Thöny-Meyer et al., 2011] Reiss, R., Ihssen, J. and Thöny-Meyer, L. (2011). Bacillus pumilus laccase: a heat stable enzyme with a wide substrate spectrum. BMC Biotechnol., 11:9.
AutoAnnotator:
Follow us:
Address:
iGEM Team TU-Munich
Emil-Erlenmeyer-Forum 5
85354 Freising, Germany
Email: igem@wzw.tum.de
Phone: +49 8161 71-4351