use of cbit.vcell.units.VCUnitDefinition in project vcell by virtualcell.
the class SpeciesContextSpec method initializeForSpatial.
public void initializeForSpatial() {
if (getDiffusionParameter() != null && getDiffusionParameter().getExpression() != null && getDiffusionParameter().getExpression().isZero()) {
Expression e = null;
ModelUnitSystem modelUnitSystem = getSimulationContext().getModel().getUnitSystem();
VCUnitDefinition micronsqpersecond = modelUnitSystem.getInstance("um2.s-1");
if (speciesContext.getStructure() instanceof Feature) {
RationalNumber rn = RationalNumber.getApproximateFraction(micronsqpersecond.convertTo(10, getDiffusionParameter().getUnitDefinition()));
e = new Expression(rn.doubleValue());
} else if (speciesContext.getStructure() instanceof Membrane) {
RationalNumber rn = RationalNumber.getApproximateFraction(micronsqpersecond.convertTo(0.1, getDiffusionParameter().getUnitDefinition()));
e = new Expression(rn.doubleValue());
} else {
RationalNumber rn = RationalNumber.getApproximateFraction(micronsqpersecond.convertTo(1.0, getDiffusionParameter().getUnitDefinition()));
e = new Expression(rn.doubleValue());
}
try {
getDiffusionParameter().setExpression(e);
} catch (ExpressionBindingException e1) {
e1.printStackTrace();
throw new RuntimeException("Error while initializing diffusion rate, " + e1.getMessage());
}
}
}
use of cbit.vcell.units.VCUnitDefinition in project vcell by virtualcell.
the class SpeciesContextSpec method getLengthPerTimeUnit.
private VCUnitDefinition getLengthPerTimeUnit() {
ModelUnitSystem modelUnitSystem = getSimulationContext().getModel().getUnitSystem();
VCUnitDefinition lengthPerTimeUnit = modelUnitSystem.getLengthUnit().divideBy(modelUnitSystem.getTimeUnit());
return lengthPerTimeUnit;
}
use of cbit.vcell.units.VCUnitDefinition in project vcell by virtualcell.
the class StochMathMapping method addJumpProcesses.
private void addJumpProcesses(VariableHash varHash, GeometryClass geometryClass, SubDomain subDomain) throws ExpressionException, ModelException, MappingException, MathException {
// set up jump processes
// get all the reactions from simulation context
// ReactionSpec[] reactionSpecs = simContext.getReactionContext().getReactionSpecs();---need to take a look here!
ModelUnitSystem modelUnitSystem = getSimulationContext().getModel().getUnitSystem();
ReactionSpec[] reactionSpecs = getSimulationContext().getReactionContext().getReactionSpecs();
for (ReactionSpec reactionSpec : reactionSpecs) {
if (reactionSpec.isExcluded()) {
continue;
}
// get the reaction
ReactionStep reactionStep = reactionSpec.getReactionStep();
Kinetics kinetics = reactionStep.getKinetics();
// probability parameter from modelUnitSystem
VCUnitDefinition probabilityParamUnit = modelUnitSystem.getStochasticSubstanceUnit().divideBy(modelUnitSystem.getTimeUnit());
// Different ways to deal with simple reactions and flux reactions
if (// simple reactions
reactionStep instanceof SimpleReaction) {
// check the reaction rate law to see if we need to decompose a reaction(reversible) into two jump processes.
// rate constants are important in calculating the probability rate.
// for Mass Action, we use KForward and KReverse,
// for General Kinetics we parse reaction rate J to see if it is in Mass Action form.
Expression forwardRate = null;
Expression reverseRate = null;
if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction) || kinetics.getKineticsDescription().equals(KineticsDescription.General)) {
Expression rateExp = new Expression(kinetics.getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate), reactionStep.getNameScope());
Parameter forwardRateParameter = null;
Parameter reverseRateParameter = null;
if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction)) {
forwardRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KForward);
reverseRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KReverse);
}
MassActionSolver.MassActionFunction maFunc = MassActionSolver.solveMassAction(forwardRateParameter, reverseRateParameter, rateExp, reactionStep);
if (maFunc.getForwardRate() == null && maFunc.getReverseRate() == null) {
throw new MappingException("Cannot generate stochastic math mapping for the reaction:" + reactionStep.getName() + "\nLooking for the rate function according to the form of k1*Reactant1^Stoir1*Reactant2^Stoir2...-k2*Product1^Stoip1*Product2^Stoip2.");
} else {
if (maFunc.getForwardRate() != null) {
forwardRate = maFunc.getForwardRate();
}
if (maFunc.getReverseRate() != null) {
reverseRate = maFunc.getReverseRate();
}
}
} else // if it's macro/microscopic kinetics, we'll have them set up as reactions with only forward rate.
if (kinetics.getKineticsDescription().equals(KineticsDescription.Macroscopic_irreversible) || kinetics.getKineticsDescription().equals(KineticsDescription.Microscopic_irreversible)) {
Expression Kon = getIdentifierSubstitutions(new Expression(reactionStep.getKinetics().getKineticsParameterFromRole(Kinetics.ROLE_KOn), getNameScope()), reactionStep.getKinetics().getKineticsParameterFromRole(Kinetics.ROLE_Binding_Radius).getUnitDefinition(), geometryClass);
if (Kon != null) {
Expression KonCopy = new Expression(Kon);
try {
MassActionSolver.substituteParameters(KonCopy, true).evaluateConstant();
forwardRate = new Expression(Kon);
} catch (ExpressionException e) {
throw new MathException(VCellErrorMessages.getMassActionSolverMessage(reactionStep.getName(), "Problem with Kon parameter in " + reactionStep.getName() + ": '" + KonCopy.infix() + "', " + e.getMessage()));
}
} else {
throw new MathException(VCellErrorMessages.getMassActionSolverMessage(reactionStep.getName(), "Kon parameter of " + reactionStep.getName() + " is null."));
}
}
boolean isForwardRatePresent = false;
boolean isReverseRatePresent = false;
if (forwardRate != null) {
isForwardRatePresent = true;
}
if (reverseRate != null) {
isReverseRatePresent = true;
}
// we process it as forward reaction
if ((isForwardRatePresent)) /*|| ((forwardRate == null) && (reverseRate == null))*/
{
// get jump process name
String jpName = TokenMangler.mangleToSName(reactionStep.getName());
// get probability
Expression exp = null;
// reactions are of mass action form
exp = getProbabilityRate(reactionStep, forwardRate, true);
ProbabilityParameter probParm = null;
try {
probParm = addProbabilityParameter(PARAMETER_PROBABILITYRATE_PREFIX + jpName, exp, PARAMETER_ROLE_P, probabilityParamUnit, reactionStep);
} catch (PropertyVetoException pve) {
pve.printStackTrace();
throw new MappingException(pve.getMessage());
}
// add probability to function or constant
varHash.addVariable(newFunctionOrConstant(getMathSymbol(probParm, geometryClass), getIdentifierSubstitutions(exp, probabilityParamUnit, geometryClass), geometryClass));
JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probParm, geometryClass)));
// actions
ReactionParticipant[] reacPart = reactionStep.getReactionParticipants();
for (int j = 0; j < reacPart.length; j++) {
Action action = null;
SpeciesCountParameter spCountParam = getSpeciesCountParameter(reacPart[j].getSpeciesContext());
if (reacPart[j] instanceof Reactant) {
// check if the reactant is a constant. If the species is a constant, there will be no action taken on this species
if (// not a constant
!simContext.getReactionContext().getSpeciesContextSpec(reacPart[j].getSpeciesContext()).isConstant()) {
int stoi = ((Reactant) reacPart[j]).getStoichiometry();
action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(-stoi));
jp.addAction(action);
}
} else if (reacPart[j] instanceof Product) {
// check if the product is a constant. If the product is a constant, there will be no action taken on this species
if (// not a constant
!simContext.getReactionContext().getSpeciesContextSpec(reacPart[j].getSpeciesContext()).isConstant()) {
int stoi = ((Product) reacPart[j]).getStoichiometry();
action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(stoi));
jp.addAction(action);
}
}
}
// add jump process to compartment subDomain
subDomain.addJumpProcess(jp);
}
if (// one more jump process for a reversible reaction
isReverseRatePresent) {
// get jump process name
String jpName = TokenMangler.mangleToSName(reactionStep.getName()) + PARAMETER_PROBABILITY_RATE_REVERSE_SUFFIX;
Expression exp = null;
// reactions are mass actions
exp = getProbabilityRate(reactionStep, reverseRate, false);
ProbabilityParameter probRevParm = null;
try {
probRevParm = addProbabilityParameter(PARAMETER_PROBABILITYRATE_PREFIX + jpName, exp, PARAMETER_ROLE_P_reverse, probabilityParamUnit, reactionStep);
} catch (PropertyVetoException pve) {
pve.printStackTrace();
throw new MappingException(pve.getMessage());
}
// add probability to function or constant
varHash.addVariable(newFunctionOrConstant(getMathSymbol(probRevParm, geometryClass), getIdentifierSubstitutions(exp, probabilityParamUnit, geometryClass), geometryClass));
JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probRevParm, geometryClass)));
// actions
ReactionParticipant[] reacPart = reactionStep.getReactionParticipants();
for (int j = 0; j < reacPart.length; j++) {
Action action = null;
SpeciesCountParameter spCountParam = getSpeciesCountParameter(reacPart[j].getSpeciesContext());
if (reacPart[j] instanceof Reactant) {
// check if the reactant is a constant. If the species is a constant, there will be no action taken on this species
if (// not a constant
!simContext.getReactionContext().getSpeciesContextSpec(reacPart[j].getSpeciesContext()).isConstant()) {
int stoi = ((Reactant) reacPart[j]).getStoichiometry();
action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(stoi));
jp.addAction(action);
}
} else if (reacPart[j] instanceof Product) {
// check if the product is a constant. If the product is a constant, there will be no action taken on this species
if (// not a constant
!simContext.getReactionContext().getSpeciesContextSpec(reacPart[j].getSpeciesContext()).isConstant()) {
int stoi = ((Product) reacPart[j]).getStoichiometry();
action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(-stoi));
jp.addAction(action);
}
}
}
// add jump process to compartment subDomain
subDomain.addJumpProcess(jp);
}
// end of if(isForwardRateNonZero), if(isReverseRateNonRate)
} else if (// flux reactions
reactionStep instanceof FluxReaction) {
// we could set jump processes for general flux rate in forms of p1*Sout + p2*Sin
if (kinetics.getKineticsDescription().equals(KineticsDescription.General) || kinetics.getKineticsDescription().equals(KineticsDescription.GeneralPermeability)) {
Expression fluxRate = new Expression(kinetics.getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate), reactionStep.getNameScope());
// we have to pass the math description para to flux solver, coz somehow math description in simulation context is not updated.
// forward and reverse rate parameters may be null
Parameter forwardRateParameter = null;
Parameter reverseRateParameter = null;
if (kinetics.getKineticsDescription().equals(KineticsDescription.GeneralPermeability)) {
forwardRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_Permeability);
reverseRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_Permeability);
}
MassActionSolver.MassActionFunction fluxFunc = MassActionSolver.solveMassAction(forwardRateParameter, reverseRateParameter, fluxRate, (FluxReaction) reactionStep);
// create jump process for forward flux if it exists.
Expression rsStructureSize = new Expression(reactionStep.getStructure().getStructureSize(), getNameScope());
VCUnitDefinition probRateUnit = modelUnitSystem.getStochasticSubstanceUnit().divideBy(modelUnitSystem.getAreaUnit()).divideBy(modelUnitSystem.getTimeUnit());
Expression rsRateUnitFactor = getUnitFactor(probRateUnit.divideBy(modelUnitSystem.getFluxReactionUnit()));
if (fluxFunc.getForwardRate() != null && !fluxFunc.getForwardRate().isZero()) {
Expression rate = fluxFunc.getForwardRate();
// get species expression (depend on structure, if mem: Species/mem_Size, if vol: species*KMOLE/vol_size)
if (fluxFunc.getReactants().size() != 1) {
throw new MappingException("Flux " + reactionStep.getName() + " should have only one reactant.");
}
SpeciesContext scReactant = fluxFunc.getReactants().get(0).getSpeciesContext();
Expression scConcExpr = new Expression(getSpeciesConcentrationParameter(scReactant), getNameScope());
Expression probExp = Expression.mult(rate, rsRateUnitFactor, rsStructureSize, scConcExpr);
// jump process name
// +"_reverse";
String jpName = TokenMangler.mangleToSName(reactionStep.getName());
ProbabilityParameter probParm = null;
try {
probParm = addProbabilityParameter(PARAMETER_PROBABILITYRATE_PREFIX + jpName, probExp, PARAMETER_ROLE_P, probabilityParamUnit, reactionStep);
} catch (PropertyVetoException pve) {
pve.printStackTrace();
throw new MappingException(pve.getMessage());
}
// add probability to function or constant
String ms = getMathSymbol(probParm, geometryClass);
Expression is = getIdentifierSubstitutions(probExp, probabilityParamUnit, geometryClass);
Variable nfoc = newFunctionOrConstant(ms, is, geometryClass);
varHash.addVariable(nfoc);
JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probParm, geometryClass)));
// actions
Action action = null;
SpeciesContext sc = fluxFunc.getReactants().get(0).getSpeciesContext();
if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(-1));
jp.addAction(action);
}
sc = fluxFunc.getProducts().get(0).getSpeciesContext();
if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(1));
jp.addAction(action);
}
subDomain.addJumpProcess(jp);
}
// create jump process for reverse flux if it exists.
if (fluxFunc.getReverseRate() != null && !fluxFunc.getReverseRate().isZero()) {
// jump process name
String jpName = TokenMangler.mangleToSName(reactionStep.getName()) + PARAMETER_PROBABILITY_RATE_REVERSE_SUFFIX;
Expression rate = fluxFunc.getReverseRate();
// get species expression (depend on structure, if mem: Species/mem_Size, if vol: species*KMOLE/vol_size)
if (fluxFunc.getProducts().size() != 1) {
throw new MappingException("Flux " + reactionStep.getName() + " should have only one product.");
}
SpeciesContext scProduct = fluxFunc.getProducts().get(0).getSpeciesContext();
Expression scConcExpr = new Expression(getSpeciesConcentrationParameter(scProduct), getNameScope());
Expression probExp = Expression.mult(rate, rsRateUnitFactor, rsStructureSize, scConcExpr);
ProbabilityParameter probRevParm = null;
try {
probRevParm = addProbabilityParameter(PARAMETER_PROBABILITYRATE_PREFIX + jpName, probExp, PARAMETER_ROLE_P_reverse, probabilityParamUnit, reactionStep);
} catch (PropertyVetoException pve) {
pve.printStackTrace();
throw new MappingException(pve.getMessage());
}
// add probability to function or constant
varHash.addVariable(newFunctionOrConstant(getMathSymbol(probRevParm, geometryClass), getIdentifierSubstitutions(probExp, probabilityParamUnit, geometryClass), geometryClass));
JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probRevParm, geometryClass)));
// actions
Action action = null;
SpeciesContext sc = fluxFunc.getReactants().get(0).getSpeciesContext();
if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(1));
jp.addAction(action);
}
sc = fluxFunc.getProducts().get(0).getSpeciesContext();
if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(-1));
jp.addAction(action);
}
subDomain.addJumpProcess(jp);
}
}
}
// end of if (simplereaction)...else if(fluxreaction)
}
// end of reaction step loop
}
use of cbit.vcell.units.VCUnitDefinition in project vcell by virtualcell.
the class StructureAnalyzer method getCorrectedRateExpression.
public Expression getCorrectedRateExpression(ReactionStep reactionStep, ReactionParticipant reactionParticipant, RateType rateType) throws Exception {
if (reactionParticipant instanceof Catalyst) {
throw new Exception("Catalyst " + reactionParticipant + " doesn't have a rate for this reaction");
// return new Expression(0.0);
}
double stoich = reactionStep.getStoichiometry(reactionParticipant.getSpeciesContext());
if (stoich == 0.0) {
return new Expression(0.0);
}
//
// make distributed rate with correct stoichiometry for this participant
//
VCUnitDefinition correctedReactionRateUnit = null;
Expression distribRate = null;
if (reactionStep.getKinetics() instanceof DistributedKinetics) {
DistributedKinetics distributedKinetics = (DistributedKinetics) reactionStep.getKinetics();
KineticsParameter distribReactionRateParameter = distributedKinetics.getReactionRateParameter();
distribRate = new Expression(distribReactionRateParameter, mathMapping.getNameScope());
correctedReactionRateUnit = distribReactionRateParameter.getUnitDefinition();
} else if (reactionStep.getKinetics() instanceof LumpedKinetics) {
//
// need to put this into concentration/time with respect to structure for reaction.
//
Structure.StructureSize structureSize = reactionStep.getStructure().getStructureSize();
LumpedKinetics lumpedKinetics = (LumpedKinetics) reactionStep.getKinetics();
KineticsParameter lumpedReactionRateParameter = lumpedKinetics.getLumpedReactionRateParameter();
Expression lumpedReactionRateExp = new Expression(lumpedReactionRateParameter, mathMapping.getNameScope());
distribRate = Expression.div(lumpedReactionRateExp, new Expression(structureSize, mathMapping.getNameScope()));
;
correctedReactionRateUnit = lumpedReactionRateParameter.getUnitDefinition().divideBy(structureSize.getUnitDefinition());
}
// correct for stoichiometry
Expression distribRateWithStoich = distribRate;
if (stoich != 1) {
distribRateWithStoich = Expression.mult(new Expression(stoich), distribRateWithStoich);
}
// flux correction if reaction and reactionParticipant are in different compartments. (not necessarily dimensionless, use KFlux parameter).
Expression distribRateWithStoichFlux = distribRateWithStoich;
if (reactionStep.getStructure() != reactionParticipant.getStructure()) {
StructureMapping reactionSM = mathMapping.getSimulationContext().getGeometryContext().getStructureMapping(reactionStep.getStructure());
StructureMapping speciesSM = mathMapping.getSimulationContext().getGeometryContext().getStructureMapping(reactionParticipant.getStructure());
Parameter fluxCorrectionParameter = mathMapping.getFluxCorrectionParameter(reactionSM, speciesSM);
Expression fluxCorrection = new Expression(fluxCorrectionParameter, mathMapping.getNameScope());
distribRateWithStoichFlux = Expression.mult(fluxCorrection, distribRateWithStoichFlux);
correctedReactionRateUnit = correctedReactionRateUnit.multiplyBy(fluxCorrectionParameter.getUnitDefinition());
}
// apply unit factor for difference substance
ModelUnitSystem unitSystem = mathMapping.getSimulationContext().getModel().getUnitSystem();
VCUnitDefinition timeUnit = unitSystem.getTimeUnit();
VCUnitDefinition speciesConcUnit = reactionParticipant.getSpeciesContext().getUnitDefinition();
VCUnitDefinition speciesConcRateUnit = speciesConcUnit.divideBy(timeUnit);
Expression unitFactor = null;
if (rateType == RateType.ConcentrationRate) {
unitFactor = mathMapping.getUnitFactor(speciesConcRateUnit.divideBy(correctedReactionRateUnit));
} else if (rateType == RateType.ResolvedFluxRate) {
unitFactor = mathMapping.getUnitFactor(speciesConcRateUnit.multiplyBy(unitSystem.getLengthUnit()).divideBy(correctedReactionRateUnit));
}
return Expression.mult(unitFactor, distribRateWithStoichFlux).flatten();
}
use of cbit.vcell.units.VCUnitDefinition in project vcell by virtualcell.
the class StructureAnalyzer method refreshTotalDependancies.
/**
* This method was created by a SmartGuide.
* @param b cbit.vcell.math.Matrix
* @param vars java.lang.String[]
*/
private void refreshTotalDependancies() throws Exception {
//
for (int v = 0; v < speciesContextMappings.length; v++) {
speciesContextMappings[v].setDependencyExpression(null);
}
if (totalNullSpaceMatrix == null) {
// System.out.println("the matrix has full rank, there are no dependencies");
return;
}
if (speciesContextMappings.length != totalNullSpaceMatrix.getNumCols()) {
throw new Exception("varName array not same dimension as b matrix");
}
StructureAnalyzer.Dependency[] dependencies = refreshTotalDependancies(totalNullSpaceMatrix, speciesContextMappings, mathMapping, false);
VCUnitDefinition totalMassUnit = null;
boolean bIsSpatial = mathMapping.getSimulationContext().getGeometry().getDimension() > 0;
ModelUnitSystem modelUnitSystem = mathMapping.getSimulationContext().getModel().getUnitSystem();
if (this instanceof VolumeStructureAnalyzer) {
if (!bIsSpatial) {
// VCUnitDefinition.UNIT_umol_um3_per_L; -> VCell vol substance unit
totalMassUnit = modelUnitSystem.getVolumeSubstanceUnit();
} else {
// VCUnitDefinition.UNIT_uM;
totalMassUnit = modelUnitSystem.getVolumeConcentrationUnit();
}
} else if (this instanceof MembraneStructureAnalyzer) {
if (!bIsSpatial) {
// VCUnitDefinition.UNIT_molecules;
totalMassUnit = modelUnitSystem.getMembraneSubstanceUnit();
} else {
// VCUnitDefinition.UNIT_molecules_per_um2;
totalMassUnit = modelUnitSystem.getMembraneConcentrationUnit();
}
}
for (int i = 0; i < dependencies.length; i++) {
String constantName = dependencies[i].invariantSymbolName;
Expression constantExp = dependencies[i].conservedMoietyExpression;
SpeciesContextMapping firstSCM = dependencies[i].speciesContextMapping;
Expression exp = dependencies[i].dependencyExpression;
//
// store totalMass parameter (e.g. K_xyz_total = xyz_init + wzy_init)
//
GeometryClass geometryClass = mathMapping.getSimulationContext().getGeometryContext().getStructureMapping(firstSCM.getSpeciesContext().getStructure()).getGeometryClass();
MathMappingParameter totalMassParameter = mathMapping.addMathMappingParameter(constantName, constantExp.flatten(), DiffEquMathMapping.PARAMETER_ROLE_TOTALMASS, totalMassUnit, geometryClass);
//
// store dependency parameter (e.g. xyz = K_xyz_total - wzy)
//
exp.bindExpression(mathMapping);
firstSCM.setDependencyExpression(exp);
}
}
Aggregations