Matlab Scripts and Functions for Chemical and Biochemical Education

This set of programs includes material created for general chemistry and biochemistry courses. They are either visual demonstrations of chemical and physical concepts or simple simulations of fundamental models in chemical and biochemical kinetics and thermodynamics. All of these programs have been tested under R2014b, Mac OS X El Capitan, but I think they should work with most versions of Matlab and most OS's. For University of Maryland students, Matlab can be downloaded via Terpware. Some of these programs make use of external functions available for download at the MathWorks¨ File Exchange, http://www.mathworks.com/matlabcentral/fileexchange/. I also provide a very simple and specialized Matlab tutorial meant to enable students to run and use canned programs, as opposed to the usual tutorials that lead them into programming Matlab. Making this material widely available was supported and inspired by the Elevate Fellows program run through the Teaching and Learning Transformation Center here at the University of Maryland College Park.

External functions export_fig, partitions, mtit, and dsxy2figxy from the Matlab File Exchange of useful user-contributed programs are needed for some of my programs. You may need to create a Matlab account to get these files. I am unsure of the protocol for providing them here so I am not doing it.

export_fig: http://www.mathworks.com/matlabcentral/fileexchange/23629-export-fig

partitions: http://www.mathworks.com/matlabcentral/fileexchange/12009-partitions-of-an-integer

mtit: http://www.mathworks.com/matlabcentral/fileexchange/3218-mtit--a-pedestrian-major-title-creator

dsxy2figxy: http://www.mathworks.com/matlabcentral/fileexchange/30347-sigplot/content/sigplot/sigplot/BasicFunctions/dsxy2figxy.m. Note this file is not findable by searching within the File Exchange -- it is within the Sigplot package. Just Google for it if necessary.

Each of my functions or scripts has a comment section (in green type) telling the user what it does and expects. They should be self-explanatory.

The Download link goes to a zip file that includes all of the files needed to run the scripts except for the external ones above. Sometimes it's just one file.

Comments welcome. If you have suggestions for improvement, or better yet improved code, please let me know.

See also my companion page of Excel spreadsheets.

Back to Jason Kahn's home page or listing of other teaching resources

Matlab programs:

Acid-Base Titration Curves. Download zip file. Back to List.

function titr = weakweak5(CA,CB,volA,volB,pKAA,pKBB,hhplot)

% Plots a titration curve for titrating a weak acid with a weak base
% Outputs a matrix with columns volume of B,pH,[HA],[A-],[B],[BH+],buffer capacity
% Jason Kahn U Maryland 4/6/2017 version 1.02 (derived from weakweak4)
% [titr = ]weakweak5(CA,CB,volA,volB,pKAA,pKBB,hhplot); for example:
% titr = weakweak5(.1,.05,100,225,4.75,3.5,0); % or just
% weakweak5(.75,1.0,300,250,6.15,-5,0); % for titration with a strong base
% if hhplot is set to 1, will also plot the pH derived from the
% Henderson-Hasselbalch just assuming complete neutralization to give [HA]
% and [A-].
% Titration of a weak acid A with concentration CA and volume volA with a
% weak base B of concentration (CB), from a volume of zero up to volB.
% KAA = 10^(-pKAA) is the acid dissociation constant of A and KBB = 10^(-pKBB)
% is the base dissociation constant of B. Volume units do not matter (as
% long as they are the same!), typically they would be ml.
% Returns a matrix with columns volB, pH,[HA],[A-],[B],[BH+],buffer capacity
% in the "titr" variable.
% The program graphs the first two versus each other to show a typical titration.
% Save the values returned in titr if you need to do multiple plots with different x
% scales. Otherwise the default is to overlay plots, so to do separate
% plots with different axes, give the "figure" command between calls to
% weakweak5.
% For example, to overlay curves showing the effect of dilution
% titr1 = weakweak5(.1,.03,100,400,4.75,3.5,1);
% titr2 = weakweak5(.01,.003,100,400,4.75,3.5,1);
% titr3 = weakweak5(.001,.0003,100,400,4.75,3.5,1);
% titr4 = weakweak5(.0001,.00003,100,400,4.75,3.5,1);
% To plot the distribution of HA, A-, B, BH+ species use these commands:
% >> titr2 = weakweak5(.1,.05,100,225,4.75,3.5,1);
% >> figure % So you don't plot on top of the pH plot
% >> plot(titr2(:,1),titr2(:,3),titr2(:,1),titr2(:,4),titr2(:,1),titr2(:,5),titr2(:,1),titr2(:,6))
% >> legend('HA','A-','B','BH+')

Acid-Base Titration Curves, slightly whizzier. Download zip file. Back to List.

function titr = weakweak7(CA,CB,volA,volB,pKAA,pKBB,hhplot,forrev,numlabel,'DualPlot')

% Plots a titration curve for titrating a weak acid with a weak base, mathematically identical to above but with more display options

Diprotic Acid Titration Curves. Download zip file. Back to List.

function titr = diprotic1(CA,CB,volA,volB,pKAA1,pKAA2,pKBB,hhplot) This version just plots the tritration.

function titr = diprotic2(CA,CB,volA,volB,pKAA1,pKAA2,pKBB,hhplot) This version plots speciation curves as well.

% Plots a titration curve for titrating a weak acid with a weak base
% Outputs a matrix with columns being:
% volume of B,pH,[H2A],[HA-],[A2-],[B],[BH+],buffer capacity
% Jason Kahn U Maryland 10/18/18 version 1.2 (derived from weakweak4)
% [titr = ]diprotic1(CA,CB,volA,volB,pKAA1,pKAA2,pKBB,hhplot); for example:
% titr = diprotic1(.1,.05,100,500,2.75,6.31,3.5,0); % or just
% diprotic1(.75,1.0,200,450,3.10,6.15,-5,0); % for titration with a strong base
% if hhplot is set to 1, will also plot the pH derived from the
% Henderson-Hasselbalch just assuming complete neutralization to give [HA-]
% and [A2-] in each segment.
% This was all independently programmed, but it is similar to
% http://www.chembuddy.com/?left=BATE&right=pH-calculator
% Titration of a weak acid A with concentration CA and volume volA with a
% weak base B of concentration (CB), from a volume of zero up to volB.
% KAA1 = 10^(-pKAA1) is the first acid dissociation constant of A,
% KAA2 = 10^(-pKAA2) is the second acid dissociation constant of A, which must
% be greater than KAA1, and KBB = 10^(-pKBB) is the base dissociation constant of B
% Volume units do not matter (as long as they are the same!), typically they would be ml.
% Returns a matrix with columns volB, pH,[H2A],[HA],[A-],[B],[BH+],bufcap in the "titr" variable.
% The program graphs the first two versus each other to show a typical titration.
% Save the values returned in titr if you need to do multiple plots with different x
% scales or look at speciation. Otherwise the default is to overlay plots, so to do separate
% plots with different axes, give the "figure" command between calls to
% diprotic1.
% For example, to overlay curves showing the effect of dilution
% titr1 = diprotic1(.1,.03,100,800,2.75,5.22,3.5,1);
% titr2 = diprotic1(.01,.003,100,800,2.75,5.22,3.5,1);
% titr3 = diprotic1(.001,.0003,100,800,2.75,5.22,3.5,1);
% titr4 = diprotic1(.0001,.00003,100,800,2.75,5.22,3.5,1);
% To plot the distribution of H2A,HA-, A2- species use these commands:
% >> titr2 = diprotic1(.1,.05,100,500,2.75,6.31,3.5,0);
% >> figure % So you don't plot on top of the pH plot
% >> plot(titr2(:,1),titr2(:,3),titr2(:,1),titr2(:,4),titr2(:,1),titr2(:,5))
% >> legend('H2A','HA-','A2-')
% (and similarly can plot B and BH+)
% To plot the distribution of H2A,HA-, A2- species as a function of pH:
% >> titr2 = diprotic1(.1,.05,100,500,2.75,6.31,3.5,0);
% >> figure % So you don't plot on top of the pH plot
% >> plot(titr2(:,2),titr2(:,3),titr2(:,2),titr2(:,4),titr2(:,2),titr2(:,5))
% >> legend('H2A','HA-','A2-')
% (and similarly can plot B and BH+)
% To plot the above distribution and buffer capacity as a double-y:
% figure
% yyaxis left
% plot(titr2(:,2),titr2(:,3),titr2(:,2),titr2(:,4),titr2(:,2),titr2(:,5))
% yyaxis right
% plot(titr2(:,2),titr2(:,8))
% legend('H2A','HA-','A2-','Bufcap')

Chemical Thermodynamics: Illustrate relationships among ΔG, ΔG°, Q, K, and molar G’s. Download zip file. Back to List.

function [ GvsX ] = TotalG(G0A,G0B,varargin)

% TotalG(G0A,G0B,totmoles,expf);, last two arguments are optional
%   Calculates free energy of a mixture of A <-> B as a function of mole
%      fraction, plots several graphs showing patterns
% Jason Kahn U Maryland 2/25/2016 version 1.01
%   Input G0A = standard state molar free energy of A, G0B for B
%      totmoles = total # of moles. Default is 1.
%   Set expf = 'y' if you want to call export_fig to output, default is no.
%   The molar free energy of each species is the standard state free
%      energy plus RT log X. Units of G are in RT for simplicity
%   X is the mole fraction of B, so 1-X is the mole fraction of A
%   Calculation assumes volume of 1 L, so 1 mole of A corresponds to 1 M A etc.
%   Returns a matrix of Total G versus X
%   Uses export_fig to make images.
%   Download export_fig from http://www.mathworks.com/matlabcentral/fileexchange/23629-exportfig
%   Calls mtit function to place overall title above plots
%   DL mtit from http://www.mathworks.com/matlabcentral/fileexchange/3218-mtit--a-pedestrian-major-title-creator
%   Also requires dsxy2figxy function, which makes arrows, from
%   https://www.mathworks.com/matlabcentral/fileexchange/30347-sigplot/content/sigplot/sigplot/BasicFunctions/dsxy2figxy.m

Chemical Kinetics I (via numerical integration). Download zip file. Back to List.

function [ T1,Y1 ] = ABC_eq_kin(varargin)

%ABC_eq_kin Simulates kinetics of the reaction scheme below via numerical
%  integration of differential rate laws
%   A <--> B <--> C
% Jason Kahn U Maryland 9/27/2016 version 1.01
% Arguments (all optional, but order is fixed and no skipping) in order are
%    A0,B0,C0,k1,k2,k3,k4,t_end,savefig
% Returns a vector of times and an array of []'s of A, B, C
% Makes one plot, a large plain-vanilla concentration vs. time plot
% Typical call with default rate constants and initial concentrations:
%    [t,Y]= ABC_eq_kin(1,0,0,6,9,2,1,3,0);
%   ABC_eq_kin(1,0,0,60,30,2,1,2.5,0); gives very fast initial
%    equilibration and slow second step
%    ABC_eq_kin(1,0,0,6,4,2,1,5,0); intermediate behavior
%    ABC_eq_kin(1,0,0,6,0,2,0,3,0); setting reverse rates = 0 gives
%    classic A->B->C behavior
%    ABC_eq_kin(1,0,0,.6,.4,200,100,5,0); shows first step rate limiting
%    and effecitvely instantaneous equilibration of B and C
%    ABC_eq_kin(1,0,0,.6,.4,1,.5,8,0); illustrates lag in appearance of C
%   Calls function ABC_eq_kin_deltas (at the end of this file) to provide differential rates
%   Calls function export_fig if and only if savefig flag is set to 1
%    Download export_fig from http://www.mathworks.com/matlabcentral/fileexchange/23629-exportfig

Chemical Kinetics II (via numerical integration). Download zip file. Back to List.

function [ T1,Y1 ] = ABCDEFkin_2(varargin)

%ABCDEFkin_2 Simulates kinetics of the reaction scheme below via numerical
%   integration of differential rate laws
%   A + B <-> C + D, D -> E, C <-> F
% Jason Kahn U Maryland 10/21/2016 version 1.1
% Arguments (all optional, but order is fixed and no skipping) in order are
%     A0,B0,C0,D0,E0,F0,k1,k2,k3,k4,k5,t_end,savefig
% Returns a vector of times and an array of []'s of A, B, C, D, E, F
% Makes two plots, one large to show all []'s and one small to evaluate SSA
%     and rapid equilibrium approximations for D and F
% Typical call with default rate constants and initial concentrations:
%     [t,Y]= ABCDEFkin_2(1,.75,0,0,0,0,6,9,2,1,.1,5,0);
% This choice of k's demos a situation where the SSA applies to D and rapid eq between C and F
%     [t,Y]= ABCDEFkin_2(1,.75,0,0,0,0,6,9,200,1000,100,2,0);
% This shows a simple A + B -> C + D -> C + E and leaves off the savefig flag
%     [t,Y]= ABCDEFkin_2(1,1,0,0,0,0,6,0,2,0,0,2);
% Calls function ABCDEFdeltas (at the end of this file) to provide differential rates
% Calls function export_fig if and only if savefig flag is set to 1
%   Download export_fig from http://www.mathworks.com/matlabcentral/fileexchange/23629-exportfig
%
% Define rate constants and initial concentrations
% Defaults:

Microstates and Configurations Animation. Download zip file. Back to List. The zip file includes several other necessary routines, but you will need to download the partitions function from the Matlab file exchange as described above.

function [ W ] = micro_mov_r3(varargin)

function outstruct = micro_movr3(varargin)
% Usage: outputstructure = micro_movr3(quanta,particles,pcflag,toppc,nodispflag,inputstructure);
% Illustrates and animates configurations and microstates for a toy model
% of distributing energy among particles, inspired by Nash's ChemThermo
% and equal to a 1-D simplified version of the Einstein model for a solid.
% Calculates W, T, S, and other results.
% Direct enumeration is SLOW for big q and p numbers! Set the nodispflag
% flag to 1 to forego looking at individual configurations and simply
% output values calculated from the Einstein solid model formulae.
% Arguments: micro_movr3(quanta,particles,pcflag,nodispflag,inputstructure,toppc);
% All arguments are optional but must be given in order if uses, some
% conflict with each other.
% Default >> micro_movr3 is equivalent to >> micro_movr3(3,5,'',8,0,'');
% Set quanta and particles and allow animations as in micro_movr3(6,8);
% set pcflag to 'pc' to show only predominant configurations, as in
% micro_movr3(13,20,'pc');
% set toppc to change the number of configurations displayed and averaged
% over in 'pc' view (default toppc = 8). Takes a numerical argument
% (e.g. 8, 15, 31, 6n-1 for n >= 6) as in
% res3 = micro_movr3(17,24,'pc',15,'',''); Toppc does nothing in
% animation mode, but if p and q are too large and the program goes
% into pc mode it will use the value of toppc. So for example
% micro_movr3(17,24,'',15,'',''); will give the same result.
% set nodispflag = 0 or omit or enter [] to allow animations if possible
% set nodispflag = 1 to just get a calculation of W, S, N0, and T as in
% micro_movr3(47,20,'',1). If nodispflag is 1, then we still try count the
% total number of configurations rather than enumerating each of them,
% but we do not generate graphs of configurations, so nodispflag = 1
% overrides 'pc' or '' in the previous argument.
% Set nodispflag = 2 to get a plot of one representative predominant
% configuration. This routine assumes the Boltzmann dsitribution is
% true and calculates a configuration that is usually very close to the
% most populated. Works for very large numbers of quanta and particles.
% If desired, could add more logic that uses the actual value of nodsipflag.
% set inputstructure to be the result structure (outstruct) from a previous
% run, to save time by reusing a partition list (qlist) under some conditions
% res1 = micro_movr3(17,20,'',0); % But note if nodispflag is set to 1 here,
% there will be no qlist generated
% res2 = micro_movr3(17,24,'',0,res1); % This will find and re-use the
% qlist from the res1 structure, if it still exists
% (RFE : the partition list generation is the slowest step. Should filter
% it and write out qlists for all smaller problems. All of the configs for
% q = 50 and p = 30 are contained in the q = 50 and p = 40 list. So for
% every q,p that the user enters could just generate ql's for all smaller
% values of p. Somehow they have to get stored for re-use, which would be
% a very large file)
% RFE implement an absolutely no plots output for use within scripts. For
% now, be careful to close plots after making them or the system will bog
% down.
% Jason Kahn U Maryland 4/3/2019 version 1.3
% outputs a structure with the following fields
% outstruct.q quanta
% outstruct.p particles
% outstruct.Wd # of microstates - direct calc
% outstruct.Sd entropy - direct calc
% outstruct.NCd # of configurations - direct calc
% outstruct.N0d direct count of particles with zero energy
% outstruct.Td T from the fit to an exponential to direct calc occupancies
% outstruct.WE # of microstates - from Einstein solid (p+q-1)!/[(p-1)!q!]
% outstruct.SE entropy - from Einstein solid model
% outstruct.NCE # of configurations from integer partition function
% outstruct.N0E # of particles with zero energy in the ensemble - from
% the area of an exponential given known TE from the thermo dE/dS
% outstruct.TE T from from the thermodynamic (1/T) = dE/dS
% outstruct.ql the qlist, the list of configurations, is provided if q < p
% for possible use in a subsequent calculation, otherwise it is 0
%
% micro_movr3 calls JDK's functions placequanta, calcdist, partitionnum, partitionnumrec,
% partitionfct, pentagonalnums, and RoundSum -- they must all be in the path as well
% Also requires "partitions" function from the Matlab File Exchange:
% This one: http://www.mathworks.com/matlabcentral/fileexchange/
% 12009-partitions-of-an-integer
% Try/catches to use cprintf if it is available.
% Shows all of the configs and microstates for "quanta" units of energy
% randomly allocated into "particles" boxes.
% Animations pause often -- hit spacebar in the plot window to continue, if you hit
% the spacebar too many times will bring the front window back and you
% may need to bring plot window to the front again to continue. You
% can always control-C out from the command window.
% Can generate a lot of plots -- get rid of them with "close all"
% Caution: We are counting microstates! Will take forever for large
% numbers of particles! You're asking for trouble if you enumerate
% when either quanta or particles is > 50. In this case the program
% may decide to take shortcuts on its own.
% For just calculating W, uses expression for the Einsten model of a solid:
% Weinsolid = factorial(particles+quanta-1)/(factorial(quanta)*factorial(particles-1));
% which will be non-infinite for particles+quanta-1 < 171, otherwise
% uses some tricks to avoid evaluating the factorial and still gives W, or just
% uses the gammaln function gammaln(n) = log(factorial(n-1)) to report
% just the entropy, proportional to ln(W).
% T is in units of the Einstein Temperature (ThetaT). T is thus 1/(the exp
% decay constant) from the Boltzmann distribution.

Coin flipping illustrates the predominant configuration. Download zip file. Back to List.

function normy = nflips2

% Plots the probability of x% of Z coin tosses coming up heads
% Uses one array of Z for discrete calculations (discz) and
% one for only a continuous approximation (zgam)
% Jason Kahn U Maryland 4/20/2016 version 1.01
% Some plots are scaled to max probability ("normalized")
% Probability(H) = (Z!)/[(H!)(T!)]
% ln(P) = ln(Z!) - ln(H!) - ln((Z-H)!)
% The gamma function generalizes the factorial Gamma(X) = (X-1)!
% Could rewrite to take these arrays as arguments...

Maxwell-Boltzmann Distribution. Download zip file. Back to List.

% Mbdist_vs_mass.m script
% Script for showing the Maxwell-Boltzmann distribution for different masses
% Jason Kahn U Maryland 2/25/2016  version 1.01
% Edit the three vectors below to control illustration.
% Clumsy; should rewrite for arbitrary number of masses, and pre-compute
% maximum velocities and kinetic energies for optimal scaling
% Calls the createMBfig function to give the distributions.
% Makes two plots, one vs v and one vs K.E.

function [v,probv,vdots,pdots,figure1] = createMBfig(molwt,Tvec,maxv,verts)

% Example: [v,probv,vdots,pdots,fig1] = createMBfig(20,100:50:300,1500,1);
% Jason Kahn U Maryland 2/25/2016  version 1.01
% All arguments are required. Plain argon at RT would be createMBfig(40,298,1000,0);
% Plots Maxwell-Boltzmann distributions for P(v) vs. v for a particle
% of mass molwt, at temperatures given in Tvec, up to maximum velocity maxv
% If verts = 1 plots vertical lines for most common, average, and rms
% velocities. Returns matrices of v's and P(v)'s, and the most common
% values, and figure handles (so calling functions can delete the plot)
% Jason Kahn, 11/10/09 - modified January 2016

LeChatelier’s Principle Illustration. Download zip file. Back to List.

function [T,Y] = Sim_twoAtoBr2(varargin)

% Usage: Sim_twoAtoBr2(kf,kr,A0,B0,deltaA,deltaB,Vratio,t_end_set,saveflag)
% as in >> Sim_twoAtoBr2(3e-2,6e-2,1,4,5,5,2,0,0); % defaults
% Function to illustrate LeChatelier's principle and dynamic equilibrium by
% simulating the gas phase dimerization reaction 2 A(g) <-> B(g)
% as a function of time, with a change in volume, change in P(A), change
% in P(B), and return to original volume to illustrate induced changes
% in the rxn quotient Q and its reversion to Kp.
% Jason Kahn U Maryland 11/8/2016 version 1.3
% Change log -- 1.3 cosmetics, checked capitalization of function names
% 1.2 (r2) calculation of relaxation times to make plots
% look good
% 1.01 first version provided
% Returns a time vector and concentration matrix, although not sure what
% one would do with the returned values.
% Integrates differential rate law, plots pressures, rates, and Q vs. K
% Arguments (all optional) must be in the following order:
% kf,kr,A0,B0,deltaA,deltaB,volume change ratio (compression for >1),
% ....set relaxation time (0 for leave alone),savefig flag
% You can leave arguments off at the end but no skipping.
% With no arguments defaults to: Sim_twoAtoBr2(3e-2,6e-2,1,4,5,5,2,0,0);
% Calls function TwoAtoB (provided separately) to calculate differential
% rates that are then integrated using ode45
% Calls function A2Brates (provided separately) to regenerate rates for
% plotting purposes
% Calls function export_fig if and only if savefig flag is set to 1
% Download export_fig from http://www.mathworks.com/matlabcentral/fileexchange/23629-exportfig
%
% Define rate constants and initial partial pressures
% Defaults:

% Differential rate law for change in concs given kf and kr for 2A <-> B
% Jason Kahn U Maryland 2/25/2016  version 1.01

function dconcs = TwoAtoB(~,concs,kf,kr)
dconcs = zeros(2,1);
dconcs(1) = -2*kf*concs(1)^2 + 2*kr*concs(2);
dconcs(2) = kf*concs(1)^2 - kr*concs(2);
end

% Function provides differential rates for 2A <-> B, called by Sim_twoAtoB
% Jason Kahn U Maryland 2/25/2016  version 1.01
function A2Brates = A2Brates(concs,kf,kr)
A2Brates(1) = kf*concs(1)^2;
A2Brates(2) = kr*concs(2);
end

Michaelis-Menten Kinetics Programs and Scripts. Download zip file. Back to List.

% This is the MMgraph_scr.m script for illustrating progressively more
%   realistic enzymatic reactions.
% Jason Kahn U Maryland 2/25/2016 version 1.01
% The top row of plots is the overall reaction,
% the bottom is initial rates. It starts with pure ideal M-M kinetics,
%   then looks at the effects of product inhibition, then reversible reactions.
% All assuming steady state MM kinetics at small [E] so E does not titrate out S or P
% This is a script rather than a function to allow fiddling with all inputs
%   as desired. Calls RevMicMen(S0,P0,KmS,KmP,Keq,VmaxS,Tf), which
%   integrates under steady-state conditions.
% Print nice output with e.g. export_fig -m3 ~/Desktop/mmkin2.png if you
% have downloaded the nice export_fig from the File Exchange.
% The figure handle is returned in fig1 so export_fig(fig1, filename);
% Note this is very kludgy and options are changed throughout below!

function [ t,concs ] = RevMicMen(S0,P0,KmS,KmP,Keq,VmaxS,Tf)

% RevMicMen Calculates time courses from M-M kinetic parameters
% Jason Kahn U Maryland 2/25/2016  version 1.01
%   Outputs concs = S and P concentrations as column vectors of one matrix,
%   vs. time as a column vector
%   Typical Usage: RevMicMen(1e-4,0,1e-5,1e-3,10,10e-5,2.5);
%   By Briggs-Haldane, VmaxR = (VmaxF*KmR)/(Keq*KmF) so it is not specified
%       intependently.
%   Program integrates vf = ([S] - [P]/Keq)*(Vmaxf/(KmS + [S] + [P]*KmS/KmP))
%   See http://www.mathworks.com/help/matlab/math/ordinary-differential-equations.html
%   Syntax is [t,y] = solver(odefun,tspan,y0,options)

function [T,Y] = MMkinetics(varargin)

% Function to simulate Michaelis-Menten enzyme kinetics as a function of time,
%   assuming irreversible conversion of S to P and no product inhibition
% Jason Kahn U Maryland 2/25/2016  version 1.01
% Plots results in 6 plots comparing steady-state and pre-steady state time scales
% Arguments (all optional) must be in the following order: k1,kr1,kcat,S0,E0,t_end,savefig flag
% You can leave arguments off at the end but no skipping.
% With no arguments defaults to: MMkinetics(5e6,6e4,300,0.1,1e-5,74.7,0);
% Returns an array of times and []'s of Efree, ES, S, and P
% Calls function MMkin to provide differential rates and to regenerate rates for plotting purposes
% Calls function export_fig if and only if savefig flag is set to 1
%   Download export_fig from http://www.mathworks.com/matlabcentral/fileexchange/23629-exportfig

function dconcs = MMkin(~,concs,k1,kr1,kcat)

% Calculate the change in concentration of Efree, ES, S, and P
% Returns the change in each concentration
% Jason Kahn U Maryland 2/25/2016  version 1.01
% These deltas are integrated in the calling program.
% Assumes initial rate conditions -- no product inhibition

function [S0vector,v0,figure1] = MM_inhib(varargin)

% Usage [S0,vs,figh] = MM_inhib(Ivector,Km,KI,KIprime,Vmax,S0vector);
%   or just MM_inhib to use defaults
% Jason Kahn U Maryland 2/25/2016  version 1.01
% All arguments are optional but must be given in order
% Graphs v0 vs. [S] for each [I], and also the Lineweaver-Burke plots
% Returns S0vector, the matrix of v0 vs. S for each [I],and the figure
% handle in case you want to export it or close it.
% Assumes initial rate conditions, so no point showing a time course.
% Assumes [E] << Km and KI and KIprime so can't titrate out I
% In the workspace, can define S0vector as for example
% >> S0 = [0.05 0.07 0.1 0.3 0.5 0.75 1 1.25 1.5 2 3.5 5 7 10];
% and Ivector as >> Ivector = [0 1e-6 3e-6 1e-5 5e-5 1e-4 ];
% Defaults to MM_inhib(Ivector,2e-5,1e-5,100,6e-3,S0);
% Change the S0vector (in units of Km) to explore different ranges of [S]
% Set defaults:

function titr = diprotic2(CA,CB,volA,volB,pKAA1,pKAA2,pKBB,hhplot)