Team:IISER TVM/drylab models


Dry Lab

Mathematical modelling

Mathematical modeling is an integral part of understanding the feasibility of our project. Here, for MOLDEMORT, we have attempted to model extremely crucial components of our project using mathematical equations to predict how our chimeric protein may work in-vitro. We have provided all the equations that have been in use, along with our understanding of the requirements of the models.

The data thus obtained can be further used to design experiments in the wet lab.

Fungal growth models

The team has identified seven fungal species through morphology and Sanger sequencing.

To check the activity of our chitinase enzymes, we need to determine the growth rate and lag phase of the species to compare it with the changed parameters when inflicted with the extracted protein.

To analyze the growth of the fungus, two different methods were used:

  1. Surface area based growth model

    When grown on a Petri plate in a PDA (Potato dextrose agar) medium, fungus tends to grow radially (given no contamination). The increasing radius and hence the increasing surface area is a parameter that can be used to check its growth rate over time.[1]

    The unidentified samples were inoculated on Petri plates and observed for 4-5 days to capture images of the growing fungus.

    Note: It was assumed that inoculum for all fungal species was of equal size.

    The images were aligned and analyzed using the ImageJ software, where the radius and the corresponding area were measured.

    Fig1: Measuring radius, area using ImageJ software

    The team has considered several growth models, which were mostly used for bacterial growth, yet can be used for fungal growth as well.[2][3]




    Where, y0= Initial log count or absorbance,y= Log count or absorbance at time t, μ= Maximum growth rate, λ= Lag time,C= Increase in log count or absorbance from y0 to ymax, β= Model coefficient.

    These three equations were used to fit the curve and find the unknown parameters. It has been observed that the Gompertz model of growth gives us the least RMSE(Root mean square error), hence proving to be the best model for our data.

    Fungal sampleGOODNESS OF FITSSER-squareadjusted R- squareRMSEμ (max growth rate)(cm^2/Hr)λ (lag time)β(Hr)

    SSE= Sum of squared estimate of errors, RMSE= Root mean square error

    Note: The surface area based growth model was conducted before the completion of fungal identification. The names have been assigned after its identification.

    Fig2: Gompertz fit of Rhizopus oryzae

    Fig3: Gompertz fit of A.niger

    Fig4: Gompertz fit of N.indicum

    Few drawbacks in this method:

    • The growth of the fungus is confined by the surface area of the Petri plate itself, and once the fungus covers the plate, it starts growing volumetrically, making it challenging to analyze the growth.
    • Most of the fungus present on our campus are sporulating. The spread of spores within the plate causes several colonies to grow together, rendering the experiment useless, as different colonies superimpose on each other, and no proper radius is visible.
    • Sometimes, the fungal species may not grow radially but in a non-uniform shape, increasing the manual labour of calculating the area.
  2. OD based growth model

    Scattering of light by a replicating population over time, and hence the changing OD is a good proxy for the growth of fungus instead of measuring its changing radius.[4]

    The fungal species were grown in a 96-well plate in a PDB (Potato dextrose broth) medium, with spores suspended in 1x PBS (Phosphate buffered saline), and observations were taken 4-5 days for an interval of 6 hours at 605nm.

    The assumption is made that the spore suspension has been uniformly distributed throughout the well.

    Fig5: Growth curve of R.oryzae

    Fig6: Growth curve of A.versicolor

    Growth rate(Absorbance/hr)Lag time(hr)
    Rhizopus oryzae0.314527.69
    Aspergillus versicolor0.170142.72

    Limitations of this method:

    • The 96 well plate was kept in an incubator-cum-shaker at 28 degrees celsius, but the heat inside the incubator caused the media in the plate to evaporate and dry up. Hence, a humidity chamber needs to be made to minimize the evaporation of media from the plate.
    • These growth models would be needed to compare the changing growth rate and lag time when the fungal species is administered with our designed chitinases.
    • A slower growth rate and higher lag time would help us to support our alternate hypothesis of developing chimeric chitinases with enhanced antifungal activity than wild-type chitinases.

IC50 of an Antifungal-Nipagin

Nipagin, being a methylparaben, is an antifungal agent primarily used in cosmetics and food preservatives.[5]

Fig7: PDB file of Nipagin visualized in Chimera

We carried out an antifungal assay using Nipagin to calculate its IC50 for the fungal species identified by the team. This same experiment would be repeated with the purified protein of our interest to calculate its IC50 for different fungal species. We wish to check the applicability of the models which calculate IC50 using Nipagin, which then would be used for our protein.

IC50 of the antifungal here would indicate how much of the chemical would be needed to inhibit the in vitro growth of the given fungal species by 50%.[6]

Lower the IC50, better the efficiency of the antifungal.

An initial trial in experimenting, with antifungal concentration starting from 1000μg/mL, with serial dilutions of ½ till the final concentration of 0.488μg/mL gave us the IC50 as follows:

  • Rhizopus oryzae: 230.46μg/mL
  • Aspergillus versicolor: 288.916μg/mL

The experiment was repeated to verify the results with different mathematical methods. Since the IC50 of Nipagin for both the fungus came in about 200-300μg/mL range, the graduation of nipagin concentration was changed to 800μg/mL to 0.5μg/mL, diluting it by 100μg/mL for the initial concentrations.

The inhibition percentage has been calculated using three different methods:[7]

  1. Method 1 Where, γ= Steepness constant
  2. Method 2 [8] Where, γ = Sigmoidal coefficient
  3. Method 3 [8] Imin= Lowest Inhibition,Imax=Highest Inhibition

Aspergillus versicolor

Fig8: IC50 of Nipagin against A.versicolor using Graphpad Prism

Fig9: IC50 of Nipagin against A.versicolor using Method_1

Fig10: IC50 of Nipagin against A.versicolor using Method_2

Fig11: IC50 of Nipagin against A.versicolor using Method_3M

Method usedHill coefficientIC50 (μg/mL)
Method 11.2785223.534
Method 21.2308241.9844
Method 31.2308241.9797
Using GraphPad-1.23241.243

Rhizopus oryzae

Fig12: IC50 of Nipagin against R.oryzae using Graphpad Prism

Fig13: IC50 of Nipagin against R.oryzae using Method_1

Fig14: IC50 of Nipagin against R.oryzae using Method_2

Fig15: IC50 of Nipagin against R.oryzae using Method_3

Method usedHill coefficientIC50 (μg/mL)
Method 12.4848231.6738
Method 22.7604294.2999
Method 32.7604294.2999
Using GraphPad-2.172255.899

Fig16: Box plot for IC50 (for different methods used)

Ligand-Receptor Interaction

Chitinase-chitin interaction can be considered as a ligand-receptor interaction, where the chitin polymer acts as a ligand that binds to the active sites of the chitinase protein, acting as receptors.

The docking of the chitin octamer to the chimeric chitinases had already been demonstrated through Autodoc results, which provides us with the binding energy.

The binding energy can be used to calculate the dissociation constants,Kd using the equation:

Here, The ligand is the chitin polymer, of which the cell wall of fungi is made of.The receptor is the active site of the chitinase protein.Formation of ligand-receptor complex:

This complex is formed by the virtue of non-covalent bonds like hydrogen bonding, electrostatic bonding, etc.

The forward rate is given by, k1[L][R] Reverse rate is given by, k2[L-R], Rate of change in the concentration of the ligand-receptor complex is given by,

Considering 1:1 ratio of ligand and receptor entity, we have,

At the steady state,

The affinity of the ligand is usually represented by the dissociation constant,

Considering the total concentration of ligand and total concentration of receptor to be constant for a small period of time, we have

Substituting [2] and [3] in equation [1], we have,

[L]=Ligand concentration at any time t, [R]=Receptor concentration at any time t,[L-R] = Ligand – Receptor complex, [LT ] = Total Ligand conc. [Chitin octamer in our case] [known] = 1uM,[RT ] = Total Receptor conc.[Chitinase in our case][known] = 1uM, k1 = Rate constant for association of complex [Determined from Kd], k2 = Rate constant for dissociation of complex [Determined from Kd],Kd =k2/k1=Dissociation constant[Determined from Autodoc protein modelling data, binding energy]

Note – This equation is for modeling interactions where 1 ligand binds to 1 receptor.

Chitinase typeLigand-Receptor complex concentration(micro Molar)
Streptomyces orientalis0.8934335

The graph has been plotted using JSim Version 2 software, by putting all the required values.

Fig17: Binding activity of the recombinant chitinases, as shown against the wild-type chitinases.

From the graph and the concentration of the Ligand-Receptor complex, it can be concluded that adding binding domains from a high-efficiency chitinase to Catalytic domains of low-efficiency chitinase increases the binding efficiency of low-efficiency chitinase, and thus the recombinant chitinase has higher binding efficiency than its parent wild type chitinase.

Note: This model depicts the binding activity of the protein to chitin octamer; itThe backward does not show the antifungal activity of the engineered protein, as only chitin polymer is taken into consideration, rather than the complex cell wall of the fungus.


We tried to estimate the amount of IPTG that would be required for the optimal protein production so as to predict and thereby help with the experimental design.

We tried to estimate the amount of IPTG that would be required for the optimal protein production so as to predict and thereby help with the experimental design.

We also predicted the max amount of protein that a single E. coli cell would be able to produce and further went on to verify our predictions with experimental results.

To get IPTG vs Protein output and Time vs Max Protein output curve, we tried modeling every metabolic process that takes place inside the cell after IPTG induction that is responsible for the production of our enzyme.

The schematic diagram of metabolic processes is attached below:

Fig18: Schematic diagram of metabolic processes after IPTG induction

The gene of our interest is expressed in the BL21 DE(3) strain of E. coli bacteria, whose expression will be under the inducible promoter. The inducer here is IPTG (Isopropyl β-d-1-thiogalactopyranoside). The T7 RNA polymerase which binds to the promoter of our gene of interest and helps in transcription is further under the control of the Lac operon. The IPTG thus have to bind with both:

  1. Lac repressor produced by LacI gene in genomic DNA ( to enable the transcription of T7 RNA polymerase) and
  2. Lac repressor produced by LacI gene in plasmid DNA so as to free the site at which T7 RNA Polymerase binds and transcribes the gene of interest.

The metabolic processes for Enzyme production are described in terms of the following Ordinary differential equations (ODEs) :[9]

For the genomic part of the E. coli ,

  1. Transcription of LacI mRNA

    Decay of LacI mRNA

  2. Translation of LacI

    Decay of LacI

  3. Formation of LacI and IPTG complex

  4. Passive diffusion of IPTG across the cell membrane of E. coli

  5. Transcription of LacY mRNA

    Decay of LacY mRNA

  6. Translation of LacY

    Decay of LacY

    Formation of IPTG-LacY complex and facilitated diffusion

  7. Formation of LacY-IPTGex complex

  8. Free LacI present (that is, not bounded to IPTG), which will prevent the transcription of T7 RNA polymerase

    From here, [IPTGin] will be replaced by [IPTGPin] ; intracellular IPTG due to facilitated diffusion by LacY produced by the plasmid of E. coli due to its copy number. The [IPTG] diffused by this virtue can bind to LacI-repressor produced by the genomic as well.
  9. Production of T7 RNA polymerase mRNA

    Decay of T7 RNA polymerase mRNA

    Binding of free LacI and mRNAT7

  10. T7 RNA Polymerase production

    Considering the plasmid part of the E. coli , the copy number of pET28a is 40.

  11. Production of LacI mRNA

  12. Production of LacY mRNA

  13. Production of LacY

  14. Free LacI

  15. Binding of RNA Polymerase

  16. Decay of mRNA Chitinase

  17. Chitinase production

    Decay of the chitinase protein

[mRNA_{LacI}]mRNA of LacI
[IPTG_{in}]Intracellular IPTG
[IPTG_{ex}]Extracellular IPTG
[mRNA_{LacY}]mRNA of LacY
[LacY] LacY
[LacY - IPTG_{ex}] LacY - IPTG_extracellular complex
[LacI - IPTG_{in}] LacI - IPTG_intracellular complex
[LacI_{F}] Free(unbounded) LacI
[mRNA_{T7 RNA POL}] mRNA of T7 RNA polymerase
[T7 RNA POL] T7 RNA polymerase
[T7 RP_{B}] T7 RNA polymerase bound to promoter region
[mRNA_{chitinase}] mRNA of Chitinase
[Chitinase] Chitinase

Note: The elements for the plasmid part are represented with the same Symbols but have an additional “P” added to them to signify plasmid.

SymbolDescription Value
KtcI,KPtcIRate constant for transcription of LacI0.23nMmin^(-1)[9]
KdmI,KPdmIDegradation constant of LacI mRNA0.462 min^(-1)[9]
Ktl,KPtlRate constant for translation of LacI15 min^(-1)[9]
K{dpI},K{PdpI} Degradation constant for LacI 0.2 min^(-1)[9]
K{pdiff},K{Ppdiff} Rate constant for passive diffusion of IPTG 0.92 min^(-1)[9]
K{fdiff},K{Pfdiff} Rate constant of facilitated diffusion 6*10^4 min^(-1)[9]
K{fY},K{PfY} Forward reaction rate constant for formation of LacY - IPTG complex 0.12nM^(-1)min^(-1)[9]
K{bY},K{PbY} Backward reaction rate constant for the formation of LacY - IPTG complex 0.1 min^(-1)[9]
K{dYI},K{PdYI} Degradation constant for LacY - IPTG complex 0.2 min^(-1)[9]
K{fI},K{PfI} Forward reaction rate constant for formation of LacY - IPTG(intracellular) complex 0.12 nM^(-1)min^(-1)[9]
K{bI},K{PbI} Backward reaction rate constant for formation of LacI - IPTG(intracellular) complex 12 min^(-1)[9]
K{tcY},K{PtcY} Rate constant for transcription of LacY 0.5nMmin^(-1)[9]
K{dmY},K{PdmY} Degradation constant of LacY mRNA 0.462 min^(-1)[9]
K{tlY},K{PtlY} Rate constant of translation of LacY 30 min^(-1)[9]
K{dY},K{PdY} Degradation constant for LacY 0.2 min^(-1)[9]
K{tcT7} Rate constant of transcription of T7 RNA Pol 0.5 nMmin^(-1)[9]
K{dmT7} Degradation constant of T7 RNA Pol mRNA 0.462 min^(-1)[9]
K{aff} Affinity constant of LacI to the operator region 0.2 min^(-1)[9]
K{tlT7} Rate constant of translation of T7 RNA Pol 6.116 min^(-1)[9]
K{dpT7} Degradation constant of T7 RNA Pol 0.462 min^(-1)[9]
K{b} Binding constant of T7 RNA to Promoter 0.22 min^(-1)[9]
K{ub} Unbinding constant 28.8 min^(-1)[9]
K{leaky} Rate constant for leaky expression 0.01 nMmin^(-1)[9]
Km Michaelis Menten constant 3*10^5 nM[9]
K1 Transcription Rate constant of Chitinase mRNA 2.3658 * 10^3 nMmin^(-1)[9]
K{dm_chitinase} Degradation Rate constant for Chitinase mRNA 0.462min^-1[9]
K{tl_chitinase} Translation Rate constant for Chitinase 10 min^-1[9]
K{dp_chitinase} Degradation constant for chitinase 0.2 min^-1[9]

Fig19: Enzyme output for different IPTG concentration

Through our model, we were able to predict that at the concentration of 0.1 - 1 mM IPTG, we can obtain optimal protein production. This prediction helped us in the design of experiments concerning the expression of the protein. It must be noted that these predictions are well in line with the standard procedures that are being used for IPTG induction.

Our other area of interest was predicting the total amount of enzyme that could be produced by a single E Coli:

Through experimental data, we got to know that 2 litres of culture produce about 20 mg of protein. Assuming there are 10^10 bacteria / L [10] of culture (general assumption) we were able to calculate the amount of protein produced by a single E. coli to be ~ 256000 nM. These experimental results are very close to our prediction which is 255853.577 nM/ cell.

Fig20: Enzyme output with Time

(Note: IPTG was induced at time t = 0)


Since our chimeric Chitinase comprises domains from two different chitinases, with both of them having a binding domain, we could, in theory, predict the activity of our chitinase by analyzing the binding of the substrate with individual domains and further coupling the activity to get a combined activity.[11]

Fig21: Binding of Substrate to Enzyme

Relative to the Total Enzyme Concentration, ET

The total velocity of this system is given by:

Where,KD1= Dissociation constant for chitinase of Serratia marcescens,KD2= Dissociation constant for chitinase of Amycolatopsis orientalis,Vmax1=Ksub>cat1[E_T] and Vmax2=Kcat2[E_T]

We here analyse the activity of bacterial Chitinase combo 2 consisting of domains from two chitinases - Serratia marcescens chitinase and Streptomyces orientalischitinase. Here instead of Streptomyces orientalis we have used proxy species of the same genus Streptomyces griseus as the literature data for Streptomyces orientalis was unavailable. The dissociation constants and Vmax for both chitinases were determined from the literature survey.[12][13]

K_{D1} 3.36*10^(-5)
K_{D2} 4.43*10^(-6)
Vmax1 0.18mmol/min
Vmax2 2.4mmol/min

The Vmax for chitinase was found to be 2.3739 mmol/min, which gives us an idea of how adding domains from high activity enzymes could in theory improve the efficiency of low activity enzymes.

Fig22: Michaelis Menten Curve

Further plotting Lineweaver Burk Plot provides us with Km and Vmax for the enzyme.

Fig23: Lineweaver Burk Plot

The predicted Vmax for chitinase was found to be 3.1756 mmol/min and Km was found to be 0.026 mM.


For the delivery of our chimeric chitinase into the human system, we are looking forward to proposing/ using nanoparticles (specifically PLGA) that are suitable due to their ability for controlled drug release, improved therapeutic efficiency, and reduction in the dose of drug administered.

Upon doing assays using our engineered chimeric protein, we have found that our protein is stable at room temperature and shows promising activity.Thus, it further strengthens our proposal to use PLGA as a delivery molecule to the human body.

PLGA - Poly(lactic-co-glycolic acid) is a synthetic copolymer made up of Lactic Acid (LA) and glycolic Acid(GA) and is extensively used for controlled drug delivery(approved by FDA). In an aqueous environment, the chain undergoes a hydrolytic attack on its ester bonds and degrades into LA and GA which are endogenous compounds normally metabolized through the Krebs cycle. By changing the LA/GA ratio and molecular weight in the polymer we could easily control the degradation time of the molecule. The degradation of the PLGA microsphere typically starts at the center of the sphere and the degradation front moves outwards.

Here we are trying to model and predict the degradation kinetics of PLGA (LA/GA ratio 50:50 and Mw 41.9 KDa in human body which helps us understanding the controlled drug release in human body)

Protein release from PLGA Microsphere(MS) is an autocatalytic polymer degradation mechanism, involving 3 prominent steps:[14]

  1. Water intruding in the device through its micro/macropores, and subsequent occurring of matrix plasticization
  2. Hydrolytically initiation of PLGA degradation at chain backbone and production of acidic by-products
  3. Produced acids further catalyze polymer degradation

Assumptions made:

  1. Spherical symmetry
  2. Invariant geometry during protein unloading

PLGA degradation governs protein solubilization and release. PLGA degradation is promoted in the MS center due to the polymer autocatalytic degradation mechanism, causing a decrease in pH and enhancing the degradation process.

Where,ρ= Non-dimensional results, RMS=Average radius of Microsphere, r= Generic radius within Microsphere

Neglecting the initial hydration phase,

Where, A= Soluble acid (carbonic acid in present blood), P= Non-degraded polymer acting as substrate, A.P= Partially degraded, insoluble polymer,k1= Rate constant for forward reaction,k-1= Rate constant for backward reaction
,k2= Rate constant for final degradation

Equation (2) can be considered as a Michaelis-Menten equation. Approximating the quasi steady state, we have,

Where,[A]0= Initial acid concentration, and [A]+[P]=[A]0

The rate of soluble product generation,


estimates the polymer-acid affinity, low KM indicates high affinity

Where,mS= Mass of soluble products, MWS= Molecular weight of soluble products, VMS= Volume of a single MS

Substituting (8) in equation (6),

Where,mP= Time-dependent mass of PLGA, MWP = Number-average molecular weight of PLGA

Substituting (10) in equation (9),

Where,mMS,0=Mass of single microspore at time 0,MWP,0= PLGA molecular weight at time 0,μP= Non-dimensional mass of residual polymer, μS = Non-dimensional mass of soluble products, ωP= Non-dimensional polymer weight normalized with respect to its initial value

Where,kdeg= Degradation constant of the polymer

Polymer loss (i.e., generation of soluble acids) takes place after an induction time, tind,

Thus, equation (15) can be modified by including a heaviside function, as follows,

To account for the continuous (non-discrete) transition of µS in time, Heaviside function has been mollified as shown in the following Equation

given the initial condition of

Assuming density of PLGA is constant during degradation, equation (14) can be written as,

Note: B0,tind, ψ are all adjustable parameters.

Fig24: Kinetics of the degradation front of PLGA MS

This model predicts that our PLGA microsphere(Mw = 41.9 KDa) will undergo full disintegration in about 60 days (2 months) in physiological conditions. The rate of protein release will be in proportion with the decay rate of the microsphere and hence the dynamics of protein released could be analyzed using this model.

  1. Aunsbjerg SD, Andersen KR, Knøchel S. Real-time monitoring of fungal inhibition and morphological changes. J Microbiol Methods. 2015;119:196-202. doi:10.1016/j.mimet.2015.10.024
  2. María-Leonor Pla, Sandra Oltra, María-Dolores Esteban, Santiago Andreu, and Alfredo Palop. (2015). Comparison of Primary Models to Predict Microbial Growth by the Plate Count and Absorbance Methods.. BioMed Research
  3. Tjørve KMC, Tjørve E (2017) The use of Gompertz models in growth analyses, and new Gompertz-model approach: An addition to the Unified-Richards family. PLoS ONE 12(6): e0178691.
  4. Langvad F. A rapid and efficient method for growth measurement of filamentous fungi. J Microbiol Methods. 1999;37(1):97-100. doi:10.1016/s0167-7012(99)00053-6
  5. National Center for Biotechnology Information (2021). PubChem Compound Summary for CID 7456, Methylparaben. Retrieved October 8, 2021 from
  6. Hoetelmans RM. "IC50 versus EC50". PK-PD relationships for anti-retroviral drugs. Amsterdam: Slotervaart Hospital. Archived from the original on 2017-05-28 – via U.S. Food and Drug Administration.
  7. Sebaugh JL. Guidelines for accurate EC50/IC50 estimation. Pharm Stat. 2011;10(2):128-134. doi:10.1002/pst.426
  8. Volpe DA, Hamed SS, Zhang LK. Use of different parameters and equations for calculation of IC₅₀ values in efflux assays: potential sources of variability in IC₅₀ determination. AAPS J. 2014;16(1):172-180. doi:10.1208/s12248-013-9554-7
  9. Stamatakis M, Mantzaris NV. Comparison of deterministic and stochastic models of the lac operon genetic network [published correction appears in Biophys J. 2009 Nov 4 ;97(9):2651]. Biophys J. 2009;96(3):887-906. doi:10.1016/j.bpj.2008.10.028
  10. Benedik, Michael. (2016). Re: What is the average number of bacteria counts per mL of broth?. Retrieved from:
  11. Zarei M, Aminzadeh S, Zolgharnein H, et al. Characterization of a chitinase with antifungal activity from a native Serratia marcescens B4A. Braz J Microbiol. 2011;42(3):1017-1029. doi:10.1590/S1517-838220110003000022
  12. M. Rabeeth, A. Anitha and Geetha Srikanth, 2011. Purification of an Antifungal Endochitinase from a Potential Biocontrol Agent Streptomyces griseus. Pakistan Journal of Biological Sciences, 14: 788-797.
  13. Netti PA, Biondi M, Frigione M. Experimental Studies and Modeling of the Degradation Process of Poly(Lactic-co-Glycolic Acid) Microspheres for Sustained Protein Release. Polymers. 2020; 12(9):2042.