文献--Bacterial coexistence driven by motility and spatial competition-模型代码

Bacterial coexistence driven by motility and spatial competition-模型代码

一、说明

• 本文主要对文献《Bacterial coexistence driven by motility and spatial competition》使用的代码模型进行注释分享。
  - [1] Gude S , Pine E , Taute K M , et al. Bacterial coexistence driven by motility and spatial competition[J]. Nature, 2020.
• 代码来源于向原作者邮件获取。

二、模型简介

• 文章中的模型采用空间扩散和Mond结合建立描述大肠杆菌在固体平面的竞争关系。
  -运行环境 matlab 2020b

三、模型建立过程

• 待续
• $$\frac{\partial p_1}{\partial t}=r^{-1}\frac{\partial }{\partial r}\left(r\frac{n}{K+n}\left(D_1\frac{\partial p_1}{\partial r}-{\chi }_1{p}_1\frac{\partial n}{\partial r}\right)\right)+\mathrm{l}\mathrm{n}2\frac{n}{K+n}{\mu }_1{p}_1$$
  $$\frac{\partial p_1}{\partial t}=r^{-1}\frac{\partial }{\partial r}\left(r\frac{n}{K+n}\left(D_1\frac{\partial p_1}{\partial r}-{\chi }_1{p}_1\frac{\partial n}{\partial r}\right)\right)+\mathrm{l}\mathrm{n}2\frac{n}{K+n}{\mu }_1{p}_1$$
  $$\frac{\partial n}{\partial t}=r^{-1}\frac{\partial }{\partial r}\left(r{D}_n\frac{\partial n}{\partial r}\right)-\mathrm{l}\mathrm{n}2\frac{n}{K+n}({\mu }_1{p}_1+{\mu }_2{p}_2)$$

四、源码

• 主函数

%% define parameters
growthRates = [0.404; 0.47]; % db/h
kMonods = [1e-3; 1e-3]; % a.u.
diffusionCoefficients = [0.02; 0.02; 1.8]; % mm^2/h
chis = [0.875; 0.0]; % mm^2/(h a.u.)
lambdas = [10; 10]; % 1/mm
plateauWidths = [1.0; 1.0]; % mm
uMaxs = [0.5e-3; 0.5e-3; 1e0]; % a.u.


xmesh = 0:0.1:17.5; % mm
tspan = 0:0.01:4*24; % hours

speciesNames = {'A', 'B', 'Nutrient'};

%% numerically solve system
solution = competition(growthRates, kMonods, diffusionCoefficients, chis, lambdas, uMaxs, plateauWidths, xmesh, tspan);

%% display final population ratio and change in population ratio
disp(strcat('Final Ratio (B/A)=', num2str((sum(solution(end,:,2).*xmesh) ./ sum(solution(end,:,1).*xmesh)))));
disp(strcat('Change Ratio =', num2str((sum(solution(end,:,2).*xmesh) ./ sum(solution(end,:,1).*xmesh)) ./ (sum(solution(1,:,2).*xmesh) ./ sum(solution(1,:,1).*xmesh)))));

%% visualize solution (A, B and nutrient)
for i=1:size(solution,3)
    
    % figure dimensions
    h_fig = figure;
    set(h_fig, 'PaperUnits','centimeters');
    set(h_fig, 'PaperPosition', [0 0 4.5 4.5]);
    
    % actual plotting
    imagesc(solution(:,:,i)');

    % axis, labels, ...
    set(gca, 'FontSize', 9);
    title(speciesNames{i});
    xlabel('Time [d]');
    ylabel('Expansion [mm]');
    set(gca, 'XTick', 0:2400:length(tspan));
    set(gca, 'XTickLabel', tspan(1):1:tspan(end)/24);
    set(gca, 'YTick', 1:50:2001);
    set(gca, 'YTickLabel', 0:5:200);
    set(gca,'YDir','normal');
    axis([1 length(tspan) 1 length(xmesh)]);
    caxis([0 4.25]);
    colormap(h_fig, 'jet');
    
    %print(h_fig, '-depsc2', strcat('TODO', speciesNames{i}, '.eps'));
end

% report amount of nutrients left over at the end of the simulation
disp(strcat('Fraction of nutrients present at the end:', num2str(calculateTotalAmount(solution(end,:,3), xmesh) / calculateTotalAmount(solution(1,:,3), xmesh)), 'a.u.'));
clear

• 模型方程

function solution = competition( growthRates, kMonods, diffusionCoefficients, chis, lambdas, uMaxs, plateauWidths, xmesh, tspan )

    % symmetry of the system
    m = 1;
    
    % numerically solve pdes
    solution = pdepe(m, @(x,t,u,DuDx)pdefun(x,t,u,DuDx,growthRates,kMonods,diffusionCoefficients,chis), @(x)icfun(x,lambdas,uMaxs,plateauWidths), @bcfun, xmesh, tspan);
end

%% Keller-Segel model and diffusion + consumption of nutrient
function [c,f,s] = pdefun(x,t,u,DuDx,growthRates,kMonods,diffusionCoefficients,chis)
    c = [1; 1; 1];
    f = [(u(3)./(kMonods(1)+u(3))).*(diffusionCoefficients(1) .* DuDx(1) - chis(1) .* u(1) .* DuDx(3)); (u(3)./(kMonods(2)+u(3))).*(diffusionCoefficients(2) .* DuDx(2) - chis(2) .* u(2) .* DuDx(3)); diffusionCoefficients(3) .* DuDx(3)];
    s = [(u(3)./(kMonods(1)+u(3))).*log(2).*growthRates(1).*u(1); (u(3)./(kMonods(2)+u(3))).*log(2).*growthRates(2).*u(2); -(log(2).*growthRates(1).*u(1).*(u(3)./(kMonods(1)+u(3))) + log(2).*growthRates(2).*u(2).*(u(3)./(kMonods(2)+u(3))))];
end

%% exponentially decaying initial profile
function u0 = icfun(x, lambdas, uMaxs, plateauWidths)
    u0 = [uMaxs(1).*exp(-lambdas(1).*(x-plateauWidths(1))); uMaxs(2).*exp(-lambdas(2).*(x-plateauWidths(2))); uMaxs(3).*ones(size(x))];
    u0(1, 1:find(x<=plateauWidths(1), 1, 'last')) = uMaxs(1);
    u0(2, 1:find(x<=plateauWidths(2), 1, 'last')) = uMaxs(2);
end

%% no flux boundary contitions on both sides
function [pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)
    pl = [0; 0; 0];
    ql = [1; 1; 1];
    pr = [0; 0; 0];
    qr = [1; 1; 1];
end

• 其他

function totalAmount = calculateTotalAmount( u , xmesh )
    
    totalAmount = u * xmesh';

end

五、文献解读

• 待续