Colloidal binding model with SingleGRM in MATLAB

Dear all,

Unfortunately I am having trouble when trying Colloidal binding with SingleGRM to model a three-protein system (sorry for still using CADET in MATLAB).

The Colloidal binding model specifies that “Salt component has to be non-binding (component 0)”. However, I am uncertain whether the SingleGRM requires salt consideration when used in conjunction with Colloidal binding model.

Any help is appreciated! Thank you in advance!

Extract of my MATLAB code:

    mGrm = SingleGRM(); 

	% Discretization
	mGrm.nComponents = 4; 
	mGrm.nCellsColumn = 50; 
	mGrm.nCellsParticle = 9; 
	mGrm.nBoundStates = ones(mGrm.nComponents, 1); % Number of bound states for each component

	% Initial conditions (equilibrated empty column)
	mGrm.initialBulk = [0.034693*1000 0.0 0 0]; % [mol / m^3], also used for the particle mobile phase 
	mGrm.initialSolid = [398 0.0 0 0]; % [mol / m^3] 
		
	% Transport
	mGrm.dispersionColumn          = 0.053077e-6; % [m^2 / s]
	mGrm.filmDiffusion             = [6.9e-6 6.9e-6 6.9e-6 6.9e-6]; % [m/s] 
	mGrm.diffusionParticle         = [7e-10 1.07e-13 7.07e-12 12.07e-12]; % [m^2 / s] 
	mGrm.diffusionParticleSurface  = [0.0 0 0 0]; % [m^2 / s]
	mGrm.interstitialVelocity      = 3.854e-4/0.8425; % [m/s]

	% Geometry
	mGrm.columnLength        = 0.185; % [m]
	mGrm.particleRadius      = 5e-5; % [m]
	mGrm.porosityColumn      = 0.37; % [-]
	mGrm.porosityParticle    = 0.75; % [-]
	
	% Adsorption
	mCol = ColloidalBinding();%
	%mCB.kineticBinding = true; % Quasi-stationary binding 
	mCol.phi                   = 1.6e8; 
	mCol.kappaExponent         = 1; 
    mCol.kappaFactor           = 1; 
	mCol.kappaConstant         = 1; 
	mCol.cordNum               = 6; 
	mCol.logKeqPhExponent                = [0 0 0]; 
	mCol.logKeqSaltPowerLawExponent      = [0.001 0.001 0.001]; 
	mCol.logKeqSaltPowerLawFactor        = [-0.1 -0.1 -0.1]; 
	mCol.logKeqSaltExpLawFactor          = [10 10 10]; 
	mCol.logKeqSaltExpLawExponentFactor  = [-10000 -10000 -10000]; 
	mCol.bppPhExponent                   = [0 0 0]; 
	mCol.bppSaltPowerLawExponent         = [0.1 0.1 0.1]; 
	mCol.bppSaltPowerLawFactor           = [1 1 1]; 
	mCol.bppSaltExpLawFactor             = [10 10 10]; 
	mCol.bppSaltExpLawExponentFactor     = [-10000 -10000 -10000]; 
	mCol.radius                          = [5.5e-9 5.5e-9 5.5e-9]; 
	mCol.kKin                            = [1 1 1]; 
	mCol.linearizationThreshold          = 1; 
	mCol.usePh                           = false; 
	mGrm.bindingModel = mCol;
	
	% Specify inlet profile

	% Reserve space: nSections x nComponents (a section can be thought of being a 
	% step in the process, see below)
	mGrm.constant       = zeros(3, mGrm.nComponents); 
	mGrm.linear         = zeros(3, mGrm.nComponents); 
	mGrm.quadratic      = zeros(3, mGrm.nComponents); 
	mGrm.cubic          = zeros(3, mGrm.nComponents); 

	% Section 1: Loading phase
	mGrm.constant(1,1)  = 0.034693*1000;  % [mol / m^3] component 1
	mGrm.constant(1,2)  = 0.219/150*10;   % [mol / m^3] component 2
    mGrm.constant(1,3)  = 0.515/150*10;   % [mol / m^3] component 3
    mGrm.constant(1,4)  = 0.152/150*10;   % [mol / m^3] component 4
    
	% Section 2: Washing phase (no protein feed)
	mGrm.constant(2,1)  = 0.034693*1000;  % [mol / m^3] component 1

	% Section 3: Elution phase (linear salt gradient with step at the beginning)
	mGrm.constant(3,1)  = 0.047544*1000;  % [mol / m^3] component 1
	mGrm.linear(3,1)  = (0.095735-0.047544)*1000/12000;  % [mol / (m^3 * s)] component 1  


	% Step 2: Create simulator and configure it

	% Construct and configure simulator
	sim = Simulator.create();
	sim.solutionTimes = linspace(0, 2092.2+1440+12000, 1001); % [s], time points at which solution is computed
	sim.sectionTimes = [0.0 2092.2 2092.2+1440 2092.2+1440+12000]; % [s] 
	sim.sectionContinuity = false(2,1);

	sim.model = mGrm;
	result = sim.run();

	solution = [result.solution.time, squeeze(result.solution.outlet{1})];
    %figure
    subplot(2,1,1);
	[hAx,hLine1,hLine2]=plotyy(solution(:, 1)-(2092.2+1440), solution(:, 3:end)*150,solution(:, 1)-(2092.2+1440), solution(:, 2));
	xlabel('Time [s]');
	ylabel('Concentration [g/L]');
    axis(hAx(1),[0,12000,0,5])
    axis(hAx(2),[0,12000,0,100])    
    hold on

Hey Leo,

welcome to the forum!

I don’t think any one of us here is still an expert on the Matlab interface (and it will get removed in a future update, so a migration to the Python interface is highly recommended) but I’ll do my best.

I think the SingleGRM as described here is just a convenience wrapper around the General Rate Model. So we can treat it like a regular GRM, which does not require a bound state for each component. So you can set the bound state for salt to 0.

Let me know if that works or doesn’t work.

1 Like

Hey Ronald,

Appreciate your quick response! With your suggestions I succeeded in running the simulation! Thank you for your kind offer of help! :smiley: