EXPERIMENTAL MODEL
For our project we wanted a way to determine the number of our sensor bacteria required in the starter culture of a given cheese in order to ensure an accurate detection of Listeria contamination. We also aimed to estimate the time taken for our sensing bacteria to respond to Listeria contamination. To accomplish these goals, we designed a series of ordinary differential equations (ODE’s) to model the rate at which each component of the agr system is up-regulated in response to AIP, and another set of ODE’s to model the response time of each step in our sensors detection system.
Each step in the regulation of AIP production and subsequent detection was defined as a “reaction” and each reaction was assigned a custom “kinetic law” to govern the rate at which that step is upregulated in response to AIP. This was done in the symbiology app of MatLab.
This model is based on equations derived from https://link.springer.com/content/pdf/10.1007%2Fs00285-009-0291-6.pdf . We assumed that equations that govern the up-regulation of the agr operon in Staphylococcus would be similar to those in Listeria.
Each reaction and its associated equation are shown below, see the end of this article for a table defining each variable and its units.
Two compartments were generated in symbiology, one for Listeria (called cell), and one for our sensing bacteria (called sensor).
For our Listeria compartment we modelled the following reactions:
1. DNA -> mRNA -> AgrA, AgrB, AgrC, AgrD
The above image shows the equations used to model transcription and translation of the agr operon. Whilst transcription of the agr components was assumed to occur at an equal rate, given that they are all components of the same operon regulated by one upstream promoter, each component of the agr system was assigned a unique rate equation to define its translation rate, and these are given below;
Highlighted here, in blue, is the equation used to describe the agr operon transcription. The equation format is as follow:
Note: the comma in the highlighted section represents an equals sign, and this format is used for all subsequent equations.
All terms are defined in an Abbreviation's table (see below).
agrB translation
agrD translation
agrA translation
agrC translation
2. Rate of AgrB integration into the membrane after AIP stimulus (TB)
3. Rate of AgrD association with transmembrane B after AIP stimulus
4.Rate of AgrC integration into the membrane after AIP stimulus
5. Rate of increase in the extracellular AIP concentration
6. AIP interaction with transmembrane AgrC
7. Rate of AgrA phosphorylation after AIP stimulus
8. Rate of interaction of Phosphorylated AgrA with the promoter for the agr operon (PII)
For the sensor, we used the following additional considerations
1. The movement of AIP through the cheese and subsequent sensor AgrC binding
Using a modified Stokes-Einstein equation, we aimed to estimate the rate of diffusion through cheese, depending on a number of important biophysical cheese properties. If we assume cheese can be represented as a porous matrix, in which the “pores” are fairly large (i.e. there is space for molecules to diffuse), then we can use Fick’s diffusion to describe mass transfer between these pores. This can be done by multiplying the diffusion coefficient (calculated using the Stokes-Einstein equation) by x/y (where x is equal to the tortuosity and y is equal to the porosity of the medium). The tortuosity of a medium (in our case cheese) can be described as the difficulty a molecule encounters when trying to move through the medium, whereas the porosity can be described as the ease a molecule encounters when moving through a medium. We, therefore, decided to assign the value of x (tortuosity) as a variable describing the total amount of fat, protein, and salt within the cheese (in other words, the dry mass of the cheese), as these are all the kind of molecules that AIP would encounter as obstacles to free diffusion, and thus they represent “difficulty”. For the value of y (porosity), we assigned a variable describing the water content in the cheese (in other words, the total weight of the cheese minus its dry mass), as we theorised AIP would encounter less difficulty moving through a cheese with a higher water content.
The resulting equation is as follows:
2. AgrA phosphorylation by our sensor AgrC upon AIP binding
3. Phosphorylated AgrA stimulated transcription of amilCP
Results and Discussion
Above is an image of the graph generated when our model is simulated in MatLab. Each line corresponds to an important component in the mechanism of agr upregulation in Listeria. For the proof-of-concept simulations random values were assigned to the parameters. However, so model was set up in such a way that we could have used an ensemble modelling approach by collecting possible values from sources like Bionumbers and creating a probability distribution for each parameter. These distributions could then have been used to conduct a Bayesian analysis for plausible predictions, including confidence intervals. We would have also liked to create a user input function, to allow a user to define the properties of their specific cheese. Most importantly, cheese types vary widely in water content, and our discussions with cheese makers had shown us that this is a very important consideration for Listeria detection. By exploring the effect of different biophysical condition, our model would thus allow us to determine an accurate prediction of the concentration of our sensor bacteria required for reliable Listeria detection in a user-specific starting culture. This again is an important feature of our project, requested by cheese makers in our stakeholder discussions with various companies (read more in our Human practices page).
Abbreviation's Table
Experimental model code
% PLEASE CHECK THAT ALL SPECIES AND PARAMETERS ARE LISTED IN THE CORRECT
% ORDER, also add in a section about initial production, not upregulation
Mobj = sbiomodel('listnot'); % Generates a new simbiology model with the name listnot
compObj = addcompartment(Mobj,'cell'); % generated a compartment named cell
compObj.CapacityUnits = 'liter'; % defines the units of the above compartment
Robj1 = addreaction(Mobj,'DNA -> DNA + M'); % Robj1 defines the reaction for the production of agr operon mRNA
Mobj.Species;
Mobj.Species(1).InitialAmount = 500; % Here we input the initial DNA concentration within the compartment 'cell'
Mobj.Species(1).InitialAmountUnits = 'molecule/milliliter';
Mobj.Species(2).InitialAmount = 0; % Here we set the initial concntration of mRNA
Mobj.Species(2).InitialAmountUnits = 'molecule/milliliter';
Kobj1 = addkineticlaw(Robj1,'AIP-mRNA2'); % adds a forward rate kinetic law (Kobj1)to define Robj1 to tell you the rate of change of mRNA (for the agr operon)
% The equation I have used is as follows, AIP-mRNA2(standing for AIP
% induced mRNA production) = (N*O*m)+(N*O*v*P)-(d*M)
Pobj1 = addparameter(Kobj1,'N');% N stands for number of listeria cells
Pobj1.Value = 20;
Pobj1.ValueUnits = 'mole/milliliter';
Pobj2 = addparameter(Kobj1,'m'); % m stands for basal rate of mRNA production
Pobj2.Value = 10;
Pobj2.ValueUnits = 'molecule/mole/second';
Pobj3 = addparameter(Kobj1,'v'); % v stands for mRNA transcription rate
Pobj3.Value = 5;
Pobj3.ValueUnits = 'molecule/mole/second';
Pobj4 = addparameter(Kobj1,'P'); % P stands for the proportion of up-regulated cells (we could maybe try to figure this out)
Pobj4.Value = 15;
Pobj4.ValueUnits = 'dimensionless';
Pobj5 = addparameter(Kobj1,'d'); % d stands for mRNA degradation rate
Pobj5.Value = 50;
Pobj5.ValueUnits = '1/second';
Pobj6 = addparameter(Kobj1,'O'); % O is the plasmid copy number (ie - the number of times the gene is going to appear in the cell)
Pobj6.Value = 1;
Pobj6.ValueUnits = 'dimensionless';
set(Kobj1,'ParameterVariableNames',{'m','v','P','d','O','N'}); % This is a list of things in the equation that are not present in the overall reaction (the parameters)
set(Kobj1,'SpeciesVariableNames',{'M'}); % This is a list of the things in the equation that do appear in the overall reaction
Robj2 = addreaction(Mobj,'M -> M + B');
Mobj.Species(3).InitialAmountUnits = 'molecule/milliliter';
Kobj2 = addkineticlaw(Robj2,'BTranslate1');
Pobj7 = addparameter(Kobj2,'K');
Pobj7.Value = 76;
Pobj7.ValueUnits = '1/second';
Pobj8 = addparameter(Kobj2,'tB');
Pobj8.Value = 49;
Pobj8.ValueUnits = '1/second';
Pobj9 = addparameter(Kobj2,'dB');
Pobj9.Value = 200000000000000000000000000000000;
Pobj9.ValueUnits = '1/second';
set(Kobj2,'ParameterVariableNames',{'K','tB','dB'});
set(Kobj2,'SpeciesVariableNames',{'M','B'});
Robj3 = addreaction(Mobj,'M -> M + D');
Mobj.Species(4).InitialAmountUnits = 'molecule/milliliter';
Kobj3 = addkineticlaw(Robj3,'DTranslate1');
Pobj10 = addparameter(Kobj3,'K');
Pobj10.Value = 76;
Pobj10.ValueUnits = '1/second';
Pobj11 = addparameter(Kobj3,'BD');
Pobj11.Value = 30;
Pobj11.ValueUnits = '1/second';
Pobj12 = addparameter(Kobj3,'dD');
Pobj12.Value = 20;
Pobj12.ValueUnits = '1/second';
set(Kobj3,'ParameterVariableNames',{'K','BD','dD'});
set(Kobj3,'SpeciesVariableNames',{'M','D'});
Robj4 = addreaction(Mobj,'M -> M + A');
Mobj.Species(5).InitialAmountUnits = 'molecule/milliliter';
Kobj4 = addkineticlaw(Robj4,'ATranslate2');
Pobj13 = addparameter(Kobj4,'K');
Pobj13.Value = 76;
Pobj13.ValueUnits = '1/second';
Pobj14 = addparameter(Kobj4,'AIPcA');
Pobj14.Value = 50;
Pobj14.ValueUnits = 'mlms';
Pobj15 = addparameter(Kobj4,'AIPCB');
Pobj15.Value = 50;
Pobj15.ValueUnits = 'molecule/milliliter';
Pobj16 = addparameter(Kobj4,'dAP');
Pobj16.Value = 20;
Pobj16.ValueUnits = '1/second';
Pobj17 = addparameter(Kobj4,'AP');
Pobj17.Value = 70;
Pobj17.ValueUnits = 'molecule/milliliter';
Pobj18 = addparameter(Kobj4,'dA');
Pobj18.Value = 9;
Pobj18.ValueUnits = '1/second';
set(Kobj4,'ParameterVariableNames',{'K','AIPcA','AIPCB','dAP','AP','dA'});
set(Kobj4,'SpeciesVariableNames',{'M','A'});
Robj5 = addreaction(Mobj,'M -> M + C');
Mobj.Species(6).InitialAmountUnits = 'molecule/milliliter';
Kobj5 = addkineticlaw(Robj5,'CTranslate3');
Pobj19 = addparameter(Kobj5,'K');
Pobj19.Value = 76;
Pobj19.ValueUnits = '1/second';
Pobj20 = addparameter(Kobj5,'tC');
Pobj20.Value = 49;
Pobj20.ValueUnits = '1/second';
Pobj21 = addparameter(Kobj5,'dC');
Pobj21.Value = 20;
Pobj21.ValueUnits = '1/second';
set(Kobj5,'ParameterVariableNames',{'K','tC','dC'});
set(Kobj5,'SpeciesVariableNames',{'M','C'});
Robj6 = addreaction(Mobj,'D + TB -> tD');
Mobj.Species(7).InitialAmount = 3000;
Mobj.Species(7).InitialAmountUnits = 'molecule/milliliter';
Mobj.Species(8).InitialAmountUnits = 'molecule/milliliter';
Kobj6 = addkineticlaw(Robj6,'AnchorD1');
Pobj22 = addparameter(Kobj6,'BD');
Pobj22.Value = 3;
Pobj22.ValueUnits = '1/second';
Pobj23 = addparameter(Kobj6,'dtD');
Pobj23.Value = 6;
Pobj23.ValueUnits = '1/second';
Pobj24 = addparameter(Kobj6,'DBAIP');
Pobj24.Value = 7;
Pobj24.ValueUnits = 'mlms';
set(Kobj6,'ParameterVariableNames',{'BD','dtD','DBAIP'});
set(Kobj6,'SpeciesVariableNames',{'D','TB','tD'});
Robj7 = addreaction(Mobj,'B -> TB');
Kobj7 = addkineticlaw(Robj7,'TransB');
Pobj25 = addparameter(Kobj7,'tB');
Pobj25.Value = 34;
Pobj25.ValueUnits = '1/second';
Pobj26 = addparameter(Kobj7,'dTB');
Pobj26.Value = 50;
Pobj26.ValueUnits = '1/second';
set(Kobj7,'ParameterVariableNames',{'tB','dTB'});
set(Kobj7,'SpeciesVariableNames',{'B','TB'});
Robj8 = addreaction(Mobj,'tD + TB -> AIP');
Mobj.Species(9).InitialAmountUnits = 'molecule/milliliter';
Kobj8 = addkineticlaw(Robj8,'AIPup');
Pobj27 = addparameter(Kobj8,'DBAIP');
Pobj27.Value = 67;
Pobj27.ValueUnits = 'mlms';
Pobj28 = addparameter(Kobj8,'AIPC');
Pobj28.Value = 20;
Pobj28.ValueUnits = 'mlms';
Pobj29 = addparameter(Kobj8,'TC');
Pobj29.Value = 30;
Pobj29.ValueUnits = 'molecule/milliliter';
Pobj30 = addparameter(Kobj8,'AIPnoC');
Pobj30.Value = 10;
Pobj30.ValueUnits = '1/second';
Pobj31 = addparameter(Kobj8,'AIPCB');
Pobj31.Value = 50;
Pobj31.ValueUnits = 'molecule/milliliter';
Pobj32 = addparameter(Kobj8,'dAIP');
Pobj32.Value = 2;
Pobj32.ValueUnits = '1/second';
set(Kobj8,'ParameterVariableNames',{'DBAIP','AIPC','TC','AIPnoC','AIPCB','dAIP'});
set(Kobj8,'SpeciesVariableNames',{'tD','TB','AIP'});
Robj9 = addreaction(Mobj,'C -> TC');
Mobj.Species(10).InitialAmountUnits = 'molecule/milliliter';
Kobj9 = addkineticlaw(Robj9,'TransC');
Pobj33 = addparameter(Kobj9,'tC');
Pobj33.Value = 60;
Pobj33.ValueUnits = '1/second';
Pobj34 = addparameter(Kobj9,'AIPC');
Pobj34.Value = 23;
Pobj34.ValueUnits = 'mlms';
Pobj35 = addparameter(Kobj9,'AIP');
Pobj35.ValueUnits ='molecule/milliliter';
Pobj36 = addparameter(Kobj9,'AIPnoC');
Pobj36.Value = 6;
Pobj36.ValueUnits = '1/second';
Pobj37 = addparameter(Kobj9,'AIPCB');
Pobj37.Value = 30;
Pobj37.ValueUnits = 'molecule/milliliter';
Pobj38 = addparameter(Kobj9,'dTC');
Pobj38.Value = 7;
Pobj38.ValueUnits = '1/second';
set(Kobj9,'ParameterVariableNames',{'tC','AIPC','AIP','AIPnoC','AIPCB','dTC'});
set(Kobj9,'SpeciesVariableNames',{'C','TC'});
Robj10 = addreaction(Mobj,'AIP + TC -> AIPCB');
Mobj.Species(11).InitialAmountUnits = 'molecule/milliliter';
Kobj10 = addkineticlaw(Robj10,'AIPboundreceptor5');
Pobj39 = addparameter(Kobj10,'AIPC');
Pobj.Value = 23;
Pobj39.ValueUnits ='mlms';
Pobj40 = addparameter(Kobj10,'AIPnoC');
Pobj40.Value = 30;
Pobj40.ValueUnits = '1/second';
Pobj41 = addparameter(Kobj10,'dAIPCB');
Pobj41.Value = 1;
Pobj41.ValueUnits = '1/second';
set(Kobj10,'ParameterVariableNames',{'AIPC','AIPnoC','dAIPCB'});
set(Kobj10,'SpeciesVariableNames',{'AIP','TC','AIPCB'});
Robj11 = addreaction(Mobj,'AIPCB + A -> AP');
Mobj.Species(12).InitialAmountUnits = 'molecule/milliliter';
Kobj11 = addkineticlaw(Robj11,'Aphosphorylation1');
Pobj42 = addparameter(Kobj11,'CAP');
Pobj42.Value = 2;
Pobj42.ValueUnits = 'mlms';
Pobj43 = addparameter(Kobj11,'noAP');
Pobj43.Value = 2;
Pobj43.ValueUnits = '1/second';
Pobj44 = addparameter(Kobj11,'dAP');
Pobj44.Value = 30;
Pobj44.ValueUnits = '1/second';
set(Kobj11,'ParameterVariableNames',{'CAP','noAP','dAP'});
set(Kobj11,'SpeciesVariableNames',{'AIPCB','A','AP'});
Robj12 = addreaction(Mobj,'AP + DNA -> M');
Kobj12 = addkineticlaw(Robj12,'ListeriaUP5');
Pobj45 = addparameter(Kobj12,'bPII');
Pobj45.Value = 5;
Pobj45.ValueUnits = 'prombinding';
Pobj46 = addparameter(Kobj12,'N');
Pobj46.Value = 20;
Pobj46.ValueUnits = 'mole/milliliter';
Pobj47 = addparameter(Kobj12,'O');
Pobj47.Value = 1;
Pobj47.ValueUnits = 'dimensionless';
Pobj48 = addparameter(Kobj12,'P');
Pobj48.Value = 0.5;
Pobj48.ValueUnits = 'dimensionless';
Pobj49 = addparameter(Kobj12,'dPII');
Pobj49.Value = 13;
Pobj49.ValueUnits = '1/second';
set(Kobj12,'ParameterVariableNames',{'bPII','N','O','P','dPII'});
set(Kobj12,'SpeciesVariableNames',{'AP','M'});
compObj2 = addcompartment(Mobj,'sensor');
compObj.CapacityUnits = 'liter';
Robj13 = addreaction(Mobj,'cell.AIP + sensor.Cs -> sensor.AIPSENSOR');
Mobj.Species(13).InitialAmountUnits = 'molecule/milliliter';
Mobj.Species(14).InitialAmountUnits = 'molecule/milliliter';
Kobj13 = addkineticlaw(Robj13,'AIPDIFFUSE4');
Pobj50 = addparameter(Kobj13,'Kb');
Pobj50.Value = 1.38064852*10-23;
Pobj50.ValueUnits = 'dimensionless'; % all parameters involved in the diffusion constant will be dimensionless but the input is not, see the equation for the units
Pobj51 = addparameter(Kobj13,'T');
Pobj51.Value = 37; % going to be variable the input must be kelvin
Pobj51.ValueUnits = 'dimensionless';
Pobj52 = addparameter(Kobj13,'vis');
Pobj52.Value = 20;
Pobj52.ValueUnits = 'dimensionless';
Pobj53 = addparameter(Kobj13,'radi');
Pobj53.Value = 2;
Pobj53.ValueUnits = 'dimensionless';
Pobj54 = addparameter(Kobj13,'fps');
Pobj54.Value = 100;
Pobj54.ValueUnits = 'dimensionless'; % need to decide a unit
Pobj55 = addparameter(Kobj13,'hho');
Pobj55.Value = 200;
Pobj55.ValueUnits = 'dimensionless';
Pobj56 = addparameter(Kobj13,'AIPCs');
Pobj56.Value = 5;
Pobj56.ValueUnits = 'mlms';
set(Kobj13,'ParameterVariableNames',{'Kb','T','vis','radi','fps','hho','AIPCs'});
set(Kobj13,'SpeciesVariableNames',{'AIP','Cs'});
Robj14 = addreaction(Mobj,'sensor.AIPSENSOR + sensor.As -> sensor.APs'); % need to add a degradation rate for this
Kobj14 = addkineticlaw(Robj14,'APCONCsense1');
Pobj57 = addparameter(Kobj14,'CAP');
Pobj57.Value = 23;
Pobj57.ValueUnits = 'mlms';
set(Kobj14,'ParameterVariableNames',{'CAP'});
set(Kobj14,'SpeciesVariableNames',{'AIPSENSOR','As'});
Robj15 = addreaction(Mobj,'sensor.APs -> sensor.ACP');
Kobj15 = addkineticlaw(Robj15,'ACPprod4');
Pobj58 = addparameter(Kobj15,'BU'); % this is binding divided by unbinding
Pobj58.Value = 30;
Pobj58.ValueUnits = '1/second';
Pobj59 = addparameter(Kobj15,'DNAs');
Pobj59.Value = 1;
Pobj59.ValueUnits = 'dimensionless';
set(Kobj15,'ParameterVariableNames',{'BU','DNAs'});
set(Kobj15,'SpeciesVariableNames',{'APs'});
csObj = getconfigset(Mobj); % Gets data set from above...I hope
optionsObj = get(csObj,'CompileOptions');
set(optionsObj,'UnitConversion',false); % Hope this is ok
set(csObj, 'MaximumNumberOfLogs', 50); % This is only here because it lets my model run even though I get negative species values (just for now as I am using random data)
unitlessObjects = sbioselect(Mobj,'*Units', '');
csObj.SolverType = 'ode15s';
csObj.StopTime = 500; % 500 second long
csObj.RuntimeOptions.StatesToLog = 'all';
[t,x,names] = sbiosimulate(Mobj);
plot(t,x);
legend(names);
xlabel('Time');
ylabel('Amount');
![]()
