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.

Hey Ronald,

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