微生物生长曲线拟合：公式详细推导与Matlab程序

微生物生长曲线拟合：公式详细推导与Matlab程序

引用

CSDN

1.

https://blog.csdn.net/qq_33980829/article/details/113623335

微生物生长曲线的拟合是微生物学研究中的一个重要环节，它可以帮助研究人员更好地理解微生物的生长规律。本文详细介绍了基于测量OD值拟合Logistic方程的具体方法，包括公式的推导过程、Matlab代码实现以及结果分析。

微生物生长曲线的拟合是微生物学研究中的一个重要环节，它可以帮助研究人员更好地理解微生物的生长规律。本文详细介绍了基于测量OD值拟合Logistic方程的具体方法，包括公式的推导过程、Matlab代码实现以及结果分析。

本文主要内容

• 基于测量的OD600进行生长曲线的拟合
• 详细的公式推导过程
• 完整的Matlab代码以及详细的代码注释
• 对公式中涉及的生物学意义进行说明

公式推导说明

• 公式基于公开的逻辑回归方程，本文主要是细化具体的推导过程
• 推导过程主要是给自己做个记录（个人需求）
• 如果对推导过程有问题，可以留言说明
• 推导过程与代码实际运行有出入

公式推导

Logistic模型的微分表达式

$$
\frac{1}{N_{i}}\frac{\partial N_{i}}{\partial t}=r({1-\tfrac{N_{i}}{K}}) \qquad (1)
$$

将右边移项到左边构成

$$
({1-\frac{N_{i}}{K}})*\frac{1}{N_{i}}{\partial N_{i}}=r\partial t \qquad (2)
$$

两边同时积分

$$
-(\frac{1}{K-N_{i}}-\frac{1}{N_{i}}){\partial N_{i}}=\int r{\partial t } \qquad (3)
$$

$$
ln \frac{K-N_{i}}{N_{i}}=-rt+c \qquad (4)
$$

两边同时取e的对数

$$
\frac{K-N_{i}}{N_{i}}=c*e^{-rt}\qquad (5)
$$

$$
\frac{K}{N_{i}}=c*e^{-rt}+1\qquad (6)
$$

两边取倒数，再两边乘K

$$
{N_{i}}= \frac{K}{c*e^{-rt}+1}\qquad (7)
$$

参数说明

• r：不同资源下的最大生长速率
• K：在特定资源下的特征生长速率
• c：微分常数

Matlab代码程序

%% 生长曲线拟合
clc
clear
close all
[~, ~, raw] = xlsread('C:\Users\huang\Desktop\生长.xlsx','Sheet1','C2:M16'); %导入文件
%% 创建输出变量
U_1= reshape([raw{:}],size(raw));
%% 清除临时变量
clearvars raw;
U_1(:,2:end)=U_1(:,2:end)/1000;
a_num=1;%标记数值
for i=4:size(U_1,2) %输入要拟合数据内容
    x=U_1(:,1);
    y=U_1(:,i);
    fx=@(b,x)(b(1)./(1+b(2).*exp(-b(3).*x)));
    %x=xlsread('data','a1:b111');输入第一组数据
    % data1
    
    % plot(x,y,'o','markerfacecolor','r')
    
    b=[1 0.6 4.3]; %初始迭代值 最大值 生长速率 （根据具体实验来设定，初始值在本方程拟合影响不大）
    for l=1:30 %拟合过程迭代
        b=lsqcurvefit(fx,b,x,y);
        b=nlinfit(x,y,fx,b);
    end
    n=length(y);
    disp('data1')
    SSy=var(y)*(n-1);
    y1=fx(b,x);
    RSS=(y-y1)'*(y-y1);
    Rsquare(1,a_num)=(SSy-RSS)/SSy;
    x1=linspace(min(x),max(x),300);
    y1=fx(b,x1);
    plot(x,y,'*')
    hold on
    b_num(a_num,:)=b;%保存每次拟合的系数信息
    a=strcat(num2str(b(1)),'./(1+',num2str(b(2)),'.*exp(-',num2str(b(3)),'.*x))'); %拟合的函数
    fplot(a,[min(x),max(x)],'linewidth',1)% 绘制拟合的函数
    %  plot(x1,y1,'-','linewidth',1) %绘制曲线
    %  plot(x1,y1,'b-','linewidth',2.5)
    hold on
    %  plot(x,y,'o')
    %   plot(x,y,'o','markerfacecolor','r')
    %   tex=texlabel(strcat(num2str(b(1)),'./(1+',num2str(b(2)),'.*exp(-',num2str(b(3)),'.*x))'));%公式表示
    %   text(5,1.1,tex) %文本位置
    hold on
    a_num=a_num+1;
    
end
xlabel('生长时间(h)','Interpreter','LaTex')
ylabel('OD$_{600}$','Interpreter','LaTex')
legend('0mg/L','拟合曲线0mg/L','50mg/L','拟合曲线50mg/L','100mg/L','拟合曲线100mg/L','150mg/L','拟合曲线150mg/L',...
    '200mg/L','拟合曲线200mg/L','250mg/L','拟合曲线250mg/L','300mg/L','拟合曲线300mg/L','350mg/L','拟合曲线350mg/L')
%  legend('Position',[0.75928385140126 0.341278249168987 0.0951822931443652 0.312833518330545])
xticks(min(x):2:max(x)+2); %x坐标区间范围以及间隔
hold off
%% 不同生长OD聚类
figure
% y_1=b_num(:,1);
% plot(b_num(:,1),b_num(:,3),'*')
%  ylim([0.3,max(b_num(:,3))])
%  yticks([0.3:0.1:max(b_num)]);
a=[1:8];
 [cidx, ctrs] =kmeans(b_num(:,1:2:3),3);
 plot(b_num(cidx==1,1),b_num(cidx==1,3),'r*', ...
             b_num(cidx==2,1),b_num(cidx==2,3),'b*',b_num(cidx==3,1),b_num(cidx==3,3),'y*');
         a=num2cell(a);
text(b_num(:,1),b_num(:,3)*1.04,a)
xlabel('k')
ylabel('r')
  

代码运行结果

生长曲线拟合

拟合的曲线聚类

数据展示

更新内容

xy=[8.531388889	58.7521
8.698333333	66.9321
8.865277778	74.1421
9.032222222	81.7941
9.199166667	89.1081
9.366111111	97.2281
9.533055556	101.3541
9.7	104.3661
9.867222222	112.7521
10.03388889	120.9121
10.20083333	128.1181
10.36777778	135.6201
10.53472222	145.9141
10.70166667	155.4581
10.86861111	165.0501
11.03555556	175.9421
11.2025	187.2941
11.36944444	201.4041
11.53638889	212.5121
11.70333333	223.4301
11.87027778	236.6141
12.0375	249.7401
12.20416667	258.9381
12.37111111	276.1021
12.53805556	293.0301
12.705	307.3221
12.87194444	324.8201
13.03888889	346.3221
13.20583333	362.6521
13.37277778	380.6641
13.53972222	402.5081
13.70666667	423.6661
13.87361111	448.1161
14.04055556	467.6561
14.2075	492.8501
14.37444444	513.2341
14.54138889	535.4541
14.70833333	554.8601
14.87527778	579.9001
15.04222222	602.2441
15.20916667	624.4821
15.37611111	645.3841
15.54305556	663.3361
15.71	685.4361
16.71166667	799.7681
16.87861111	815.7301
17.04555556	828.9801
17.2125	845.3301
17.37944444	855.9781
17.54638889	867.4361
17.71333333	877.0741
17.88027778	886.7281
18.04722222	895.7761
18.21416667	909.2841
18.38111111	919.7761
18.54805556	930.9421
18.715	940.8961
18.88194444	951.1941
20.05055556	1021.0361
20.2175	1042.6761
20.38444444	1060.8581
20.55138889	1073.5341
20.71833333	1070.9521
20.88527778	1059.4561
21.05222222	1048.9141
21.21916667	1035.6741
21.38611111	1026.7041
21.55305556	1017.1421
21.72	1012.1061
21.88694444	1003.3861
22.05388889	995.2201
22.22111111	986.8761
22.38777778	981.9881
22.55472222	977.9221
22.72166667	974.0181
22.88861111	973.4541
23.05555556	971.3061
23.2225	968.5741
23.38944444	968.6221
];
a_ppp=1;
%% 生长曲线拟合
% clc
% clear
% close all
% [~, ~, raw] = xlsread('B-Hem.xlsx','Channel5','B50:C200'); %导入文件
% %% 创建输出变量
% U_1= reshape([raw{:}],size(raw));
% %% 清除临时变量
% clearvars raw;
% U_1(:,2:end)=U_1(:,2:end);
yy_a=1;%误差梯度梯度下降
% [xData, yData] = prepareCurveData(U_1(:,1), U_1(:,2:end));
[xData, yData] = prepareCurveData( xy(:,1), xy(:,2 ));
for i=1:100
% Set up fittype and options.
ft = fittype( 'a./(1+b.*exp(-c.*x))', 'independent', 'x', 'dependent', 'y' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.StartPoint = [1000 0.764781737089429 0.335121584564207];
% Fit model to data.
[fitresult, gof] = fit(xData, yData, ft, opts );
dyy=fitresult(xData);
a_num=find(abs(yData-dyy)>dyy*0.003*yy_a);
if gof.rsquare<0.999 %拟合度大于0.95跳出循环
        xData(a_num)=[];
    yData(a_num)=[];
    yy_a=yy_a-0.0001;%每次降低部分误差  
    a_ppp=1+a_ppp;
else %删除偏差值
  break
end
end
%生长曲线拟合
a=0:0.1:30;%范围
dyy=fitresult(a);
dxx=a;
plot(dxx,dyy)
hold on 
plot(xy(:,1), xy(:,2 ),'o')
%比生长曲线拟合
a=0:0.1:30;%范围
dy=diff(fitresult(a),1)./fitresult(a(2:end));%差分(比生长速率)
dy1=diff(fitresult(a),1);
dy2=diff(fitresult(a),2);
dx=a;
%生长速率曲线
figure
plot(dx(1:end-1),dy1)%(生长速率)
%比生长曲线
figure
plot(dx(1:end-1),dy)
%二阶生长速率曲线
figure
plot(dx(1:end-2),dy2)%(生长速率)
max_dy=max(dy);%最大比生长速率
max_dy1=max(dy1);%最大生长速率
max_dy2=max(dy2);%最大二阶生长速率
disp('最大比生长速率为：')
disp(max_dy)
disp('最大生长速率为：')
disp(max_dy1)
disp('二阶最大生长速率为：')
disp(max_dy2)
  

生长曲线

生长速率曲线

比生长数量曲线