Team:Linkoping Sweden - 2021.igem.org
Matlab Code

-->modeling_clean
% =============================================================
% experimental data
    data.growthtime1 = (1:24); % x1 (days)

    data.BG11 = [0.0846	0.1063	0.1166	0.1897	0.2359	0.3094	0.3755	0.4504	0.5498	0.6260	0.7598	0.8935	1.0279	1.1395	1.1500	1.3256	1.4018	1.4990	1.6067	1.7638	1.8523	1.9469	2.1148	2.2579]; % y1 (OD)
    data.BRA = [0.0619	0.0898	0.0959	0.1642	0.2060	0.2450	0.3297	0.3686	0.4768	0.5080	0.5997	0.6913	0.7905	0.8672	0.9033	0.9456	1.0293	1.1339	1.2207	1.3292	1.3621	1.4325	1.5932	1.6459]; % y2 (OD)
    data.SEA = [0.0612	0.0794	0.1055	0.1630	0.1968	0.2364	0.2891	0.3477	0.4352	0.4605	0.5531	0.6456	0.7368	0.7985	0.8356	0.9609	1.0187	1.0638	1.1256	1.2153	1.2797	1.3206	1.3999	1.4945]; % y2 (OD)
    data.PHACHT = [0.0843	0.1212	0.1688	0.1688	0.2846	0.3372	0.4144	0.4827	0.5745	0.6529	0.7621	0.8712	1.0137	1.1093	1.1520	1.3014	1.3759	1.4595	1.6173	1.7094	1.8920	2.0768	2.4100	2.5833]; % y2 (OD)
    data.SIEBEN = [0.0720	0.1284	0.1155	0.1793	0.2398	0.3046	0.3584	0.4827	0.5274	0.5930	0.6953	0.7976	0.9148	0.9966	1.0533	1.1650	1.2432	1.3115	1.3619	1.5116	1.6012	1.6657	1.8095	1.8731]; % y2 (OD)
    data.SECHS = [0.0722	0.0935	0.1268	0.1366	0.1692	0.2012	0.2245	0.2594	0.3129	0.3452	0.4140	0.4827	0.5642	0.6170	0.6912	0.8075	0.8651	1.0637	1.0220	1.1448	1.2423	1.3122	1.4631	1.5332]; % y2 (OD)
% =============================================================
% variables

     e = 2.718;  % eulerische zahl
     A= 83.3;  % initial growthrate, pre-calculated in excel so that it fits to experimental start values
     C= 2.53;  % maximum amount of Bacteria (OD), pre-calculated in excel

     pHmin= 5;   % from paper, value for our Cyanobacteria
     pHmax= 11;  % from paper, value for our Cyanobacteria
     pHopt= 7.5; % from paper, value for our Cyanobacteria
     pH1 = 8;    % from our experiment
     pH2 = 7;    % from our experiment
     pH3 = 6;    % from our experiment

% =============================================================

% plot experimental data

    figure (1); clf
    plot(data.growthtime1, data.BG11,'r',data.growthtime1, data.BRA,'r',data.growthtime1, data.SEA,'r',data.growthtime1, data.PHACHT,'r',data.growthtime1, data.SIEBEN,'r',data.growthtime1,  data.SECHS,'r');
    xlim([0 30]); xlabel('x [days]'); ylabel('y [OD544]');

 % =============================================================

%Calculate  influence pH

    % RpH= growth rate under different pHs (value from day 5)

        data.pH1 = ( data.PHACHT ./ data.growthtime1 );
        RpH1 = data.pH1;
        RpH1max = max (data.pH1);

        % pH1 = 8;

        data.pH2 = (data.SIEBEN  ./ data.growthtime1 );
        RpH2 = data.pH2;
        RpH2max = max (data.pH2);
        % pH2 = 7;

        data.pH3 = (data.SECHS ./ data.growthtime1);
        RpH3 = data.pH3;
        RpH3max = max (data.pH3);
        % pH3 = 6;

       RpHmean= (RpH1max+RpH2max+RpH3max)/3;

     % RpH= c*(pH-pHmin)*(pH-pHmax))/(d*(pH-pHmin)*(pH-pHmax)-(e*((pH-pHopt)^2)))
     % solve pH equation for c,d and e

        syms c d g


            eqn1= (c*(pH1-pHmin)*(pH1-pHmax))/(d*((pH1-pHmin)*(pH1-pHmax)-g*(pH1-pHopt)^2))==0.207;
            eqn2= (c*(pH2-pHmin)*(pH2-pHmax))/(d*((pH2-pHmin)*(pH2-pHmax)-g*(pH2-pHopt)^2))==0.207;
            eqn3= (c*(pH3-pHmin)*(pH3-pHmax))/(d*((pH3-pHmin)*(pH3-pHmax)-g*(pH3-pHopt)^2))==0.207;

           sol = solve([eqn1, eqn2, eqn3], [c, d, g]);

           c= sol.c;

           d= sol.d;

           g= sol.g;

      % calculate RpH, growthrate dependent on pH

            pH= (5:11);
            lenpH= length(pH);

            RpHresults=zeros(1,7); %size depends on how many different pH we are running

             for z= 1:lenpH
                RpH=(c*(pH(z) - pHmin) * (pH(z)-pHmax))/((d * (pH(z)-pHmin) * (pH(z)-pHmax)) - (g*(pH(z)-pHopt)*(pH(z)-pHopt)));
                RpHresults(z)= RpH;


             end

        % calculate fph, growthcurve under different pH

            lenRpHresults= length(RpHresults);
            fpHresults=zeros ([121 7]);

            for w = 1:lenRpHresults
                 for t = 1:121
                    fpH= C/(1+(A*(e^(-RpHresults(w)*(t-1)))));
                    fpHresults(t,w)= fpH;
                 end
            end

        % plot results pH

        figure (2); clf
             plot(fpHresults);
             xlim([0 30]); xlabel('x [days]'); ylabel('y [OD544]');


% =============================================================

% Calculate  influence salinity

    % Rsal= growth rate under different salinity levels

     % RpH= growth rate under different pHs

        data.sal1 = (( data.BG11./ data.growthtime1));
        Rsal1 = max (data.sal1);
        sal1 = 0.027;  % 1.62 g/l /58.44 g/mol

        data.sal2 = ((data.BRA ./ data.growthtime1));
        Rsal2 = max (data.sal2);
        sal2 = 0.2566;  % 15 g/l /58.44 g/mol

        data.sal3 = ((data.SEA ./ data.growthtime1));
        Rsal3 = max (data.sal3);
        sal3 = 0.598;  % 35 g/l /58.44 g/mol

  % Rsal= (h* (sal^2))+ (i*sal)+ j
  % solve salinity equation for h,i and j

        syms h i j
        eqn4= (h* (sal1^2))+ (i*sal1)+ j==0.244;
        eqn5= (h* (sal2^2))+ (i*sal2)+ j==0.191;
        eqn6= (h* (sal3^2))+ (i*sal3)+ j==0.182;


        sol = solve([eqn4, eqn5, eqn6], [h, i, j]);
        h = sol.h;
        i = sol.i;
        j = sol.j;


     % calculate Rsal growthrate dependent on salinity

        sal= [0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1];
        lensal= length(sal);

        Rsalresults=zeros(1,11); %size depends on how many different salinities we are running

        for v= 1:lensal
           Rsal= (h* (sal(v)^2))+ (i*sal(v))+ j;
            Rsalresults(v)= Rsal;

        end

        % calculate fsal, growthcurve under different salinities

        lenRsalresults= length(Rsalresults);
        fsalresults=zeros ([121 11]);

        for u = 1:lenRsalresults
            for t = 1:121
                 fsal= C/(1+(A*(e^((-Rsalresults(u)*(t-1))))));
                 fsalresults(t,u)= fsal;
            end
        end

        figure (3); clf
             plot(fsalresults);
             xlim([0 120]); xlabel('x [days]'); ylabel('y [OD544]');

       fsalresults=zeros ([121 1]);
       Rsal= (h* (0.027^2))+ (i*0.027)+ j;
       for t = 1:121
                 fsal= C/(1+((A)*(e^((-Rsal*(t-1))))));
                 fsalresults(t)= fsal;
            end

% =============================================================

% Calculate Influence factors

    % ftotal= (o*fpHresults) + (p*fsalresults);


  % solve salinity equation for o and p
        ftotalA = 0.8067; % average day 17 Sea

        RpHresultsA=(c*(7.5 - pHmin) * (7.5-pHmax))/((d * (7.5-pHmin) * (7.5-pHmax)) - (g*(7.5-pHopt)*(7.5-pHopt)));
        fpHresultsA = C/(1+(A*(e^((-RpHresultsA *17)))));

        RsalresultsA= (h* (0.598^2))+ (i*0.598)+ j;
        fsalresultsA = C/(1+(A*(e^((-RsalresultsA *17)))));


        ftotalB = 0.672; % average day 17 pH6

        RpHresultsB=(c*(6 - pHmin) * (6-pHmax))/((d * (6-pHmin) * (6-pHmax)) - (g*(6-pHopt)*(6-pHopt)));
        fpHresultsB = C/(1+(A*(e^((-RpHresultsB)*17))));

        RsalresultsB= (h* (0.027^2))+ (i*0.027)+ j;
        fsalresultsB = C/(1+(A*(e^((-RsalresultsB*17)))));

        syms o p
        eqn7= ftotalA == (o*fpHresultsA) + (p*fsalresultsA); % pH = 7.5 salinity = 0.598
        eqn8= ftotalB == (o*fpHresultsB) + (p*fsalresultsB); % pH = 6 salinity = 0.027



        sol = solve([eqn7, eqn8], [o, p]);
        o = sol.o;
        p = sol.p;

     %solve equation for every condition we want

        pHmodel =  7; % here we can put the pH we want to model
        salmodel = 0.5; % here we can put the salinity we want to model

        RpHresultsmodel=(c*(pHmodel - pHmin) * (pHmodel-pHmax))/((d * (pHmodel-pHmin) * (pHmodel-pHmax)) - (e*(pHmodel-pHopt)*(pHmodel-pHopt)));

        Rsalresultsmodel= (h* (salmodel^2))+ (i*salmodel)+ j;

        ftotalresultsmodel=zeros ([121 1]);

        for t= 1:121
            fsalresultsmodel(t) = C/(1+(A*(e^((-Rsalresultsmodel*(t-1))))));
            fpHresultsmodel(t) = C/(1+(A*(e^((-RpHresultsmodel)*(t-1)))));
            ftotal= (o*fpHresultsmodel(t) + (p*fsalresultsmodel((t))));
            ftotalresultsmodel(t)= ftotal;
        end


 % plot results model

        figure (4); clf
             plot(ftotalresultsmodel);
             xlim([0 120]); xlabel('x [days]'); ylabel('y [OD544]');
Template:Linkoping Sweden/ModelingCode

Matlab Code
