use of cbit.vcell.model.ModelUnitSystem in project vcell by virtualcell.
the class ElectricalStimulus method setParameterValue.
/**
* This method was created by a SmartGuide.
* @param expressionString java.lang.String
* @exception java.lang.Exception The exception description.
*/
public void setParameterValue(LocalParameter parm, Expression exp) throws ExpressionException, PropertyVetoException {
Parameter p = parameterContext.getLocalParameterFromName(parm.getName());
if (p != parm) {
throw new RuntimeException("parameter " + parm.getName() + " not found");
}
Expression oldExpression = parm.getExpression();
boolean bBound = false;
try {
LocalParameter[] newLocalParameters = (LocalParameter[]) parameterContext.getLocalParameters().clone();
String[] symbols = exp.getSymbols();
Vector<String> symbolsToAdd = new Vector<String>();
for (int i = 0; symbols != null && i < symbols.length; i++) {
if (parameterContext.getEntry(symbols[i]) == null) {
symbolsToAdd.add(symbols[i]);
}
}
ModelUnitSystem modelUnitSystem = simulationContext.getModel().getUnitSystem();
for (int i = 0; i < symbolsToAdd.size(); i++) {
newLocalParameters = (LocalParameter[]) BeanUtils.addElement(newLocalParameters, parameterContext.new LocalParameter(symbolsToAdd.elementAt(i), new Expression(0.0), ElectricalStimulusParameterType.UserDefined, modelUnitSystem.getInstance_TBD(), ElectricalStimulusParameterType.UserDefined.databaseRoleTag));
}
parameterContext.setLocalParameters(newLocalParameters);
exp.bindExpression(parameterContext);
parm.setExpression(exp);
bBound = true;
} finally {
try {
if (!bBound) {
parm.setExpression(oldExpression);
}
parameterContext.cleanupParameters();
} catch (PropertyVetoException e) {
e.printStackTrace(System.out);
throw new RuntimeException(e.getMessage());
}
}
}
use of cbit.vcell.model.ModelUnitSystem in project vcell by virtualcell.
the class ElectricalStimulus method parameterVCMLSet.
/**
* This method was created by a SmartGuide.
* @param tokens java.util.StringTokenizer
* @exception java.lang.Exception The exception description.
*/
public final void parameterVCMLSet(CommentStringTokenizer tokens) throws ExpressionException, PropertyVetoException {
if (tokens == null) {
return;
}
Vector<LocalParameter> esParametersV = new Vector<LocalParameter>();
if (!tokens.nextToken().equalsIgnoreCase(VCMODL.ElectricalStimulus) || !tokens.nextToken().equalsIgnoreCase(GENERAL_PROTOCOL) || !tokens.nextToken().equalsIgnoreCase(VCMODL.BeginBlock)) {
throw new RuntimeException(ElectricalStimulus.class.getName() + ".parameterVCMLRead, unexpected token ");
}
String token = null;
ModelUnitSystem modelUnitSystem = simulationContext.getModel().getUnitSystem();
while (tokens.hasMoreTokens()) {
token = tokens.nextToken();
if (token.equalsIgnoreCase(VCMODL.EndBlock)) {
break;
}
if (token.equalsIgnoreCase(VCMODL.Parameter)) {
String roleName = tokens.nextToken();
ElectricalStimulusParameterType parameterType = null;
for (ElectricalStimulusParameterType role : ElectricalStimulusParameterType.values()) {
if (roleName.equals(role.databaseRoleTag)) {
parameterType = role;
break;
}
}
if (parameterType == null) {
throw new RuntimeException(ElectricalStimulus.class.getName() + ".parameterVCMLRead, unexpected token for roleName " + roleName);
}
String parameterName = tokens.nextToken();
Expression exp = MathFunctionDefinitions.fixFunctionSyntax(tokens);
String unitsString = tokens.nextToken();
VCUnitDefinition unitDef = modelUnitSystem.getInstance_TBD();
if (unitsString.startsWith("[")) {
while (!unitsString.endsWith("]")) {
String tempToken = tokens.nextToken();
unitsString = unitsString + " " + tempToken;
}
//
// now string starts with '[' and ends with ']'
//
unitDef = modelUnitSystem.getInstance(unitsString.substring(1, unitsString.length() - 1));
} else {
tokens.pushToken(unitsString);
}
LocalParameter esp = parameterContext.new LocalParameter(parameterName, exp, parameterType, unitDef, parameterType.databaseRoleTag);
esParametersV.add(esp);
} else {
throw new RuntimeException(ElectricalStimulus.class.getName() + ".parameterVCMLRead, unexpected token for paramter tag " + token);
}
}
if (esParametersV.size() > 0) {
LocalParameter[] espArr = new LocalParameter[esParametersV.size()];
esParametersV.copyInto(espArr);
parameterContext.setLocalParameters(espArr);
} else {
parameterContext.setLocalParameters(null);
}
}
use of cbit.vcell.model.ModelUnitSystem in project vcell by virtualcell.
the class MembraneStructureAnalyzer method refreshResolvedFluxes.
/**
* This method was created in VisualAge.
*/
private void refreshResolvedFluxes() throws Exception {
// System.out.println("MembraneStructureAnalyzer.refreshResolvedFluxes()");
ModelUnitSystem unitSystem = mathMapping.getSimulationContext().getModel().getUnitSystem();
GeometryContext geoContext = mathMapping.getSimulationContext().getGeometryContext();
Vector<ResolvedFlux> resolvedFluxList = new Vector<ResolvedFlux>();
//
// for each reaction, get all fluxReactions associated with this membrane
//
Vector<FluxReaction> fluxList = new Vector<FluxReaction>();
ReactionSpec[] reactionSpecs = mathMapping.getSimulationContext().getReactionContext().getReactionSpecs();
for (int j = 0; j < reactionSpecs.length; j++) {
if (reactionSpecs[j].isExcluded()) {
continue;
}
ReactionStep rs = reactionSpecs[j].getReactionStep();
if (rs.getStructure() != null && geoContext.getStructureMapping(rs.getStructure()).getGeometryClass() == surfaceClass) {
if (rs instanceof FluxReaction) {
fluxList.addElement((FluxReaction) rs);
}
}
}
//
for (int i = 0; i < fluxList.size(); i++) {
FluxReaction fr = fluxList.elementAt(i);
ReactionParticipant[] reactionParticipants = fr.getReactionParticipants();
for (int j = 0; j < reactionParticipants.length; j++) {
if (!(reactionParticipants[j] instanceof Reactant) && !(reactionParticipants[j] instanceof Product)) {
continue;
}
ResolvedFlux rf = null;
SpeciesContext speciesContext = reactionParticipants[j].getSpeciesContext();
for (int k = 0; k < resolvedFluxList.size(); k++) {
ResolvedFlux rf_tmp = resolvedFluxList.elementAt(k);
if (rf_tmp.getSpeciesContext() == reactionParticipants[j].getSpeciesContext()) {
rf = rf_tmp;
}
}
//
// if speciesContext is not "fixed" and is mapped to a volume, add flux to ResolvedFlux
//
StructureMapping structureMapping = mathMapping.getSimulationContext().getGeometryContext().getStructureMapping(reactionParticipants[j].getStructure());
if (structureMapping.getGeometryClass() == surfaceClass) {
// flux within surface
continue;
}
if (structureMapping.getGeometryClass() instanceof SubVolume && surfaceClass.isAdjacentTo((SubVolume) structureMapping.getGeometryClass())) {
SpeciesContextSpec speciesContextSpec = mathMapping.getSimulationContext().getReactionContext().getSpeciesContextSpec(speciesContext);
if (!speciesContextSpec.isConstant()) {
if (rf == null) {
VCUnitDefinition speciesFluxUnit = speciesContext.getUnitDefinition().multiplyBy(unitSystem.getLengthUnit()).divideBy(unitSystem.getTimeUnit());
rf = new ResolvedFlux(speciesContext, speciesFluxUnit);
resolvedFluxList.addElement(rf);
}
FeatureMapping featureMapping = (FeatureMapping) structureMapping;
Expression insideFluxCorrection = Expression.invert(new Expression(featureMapping.getVolumePerUnitVolumeParameter(), mathMapping.getNameScope()));
//
if (fr.getKinetics() instanceof DistributedKinetics) {
KineticsParameter reactionRateParameter = ((DistributedKinetics) fr.getKinetics()).getReactionRateParameter();
Expression correctedReactionRate = Expression.mult(new Expression(reactionRateParameter, mathMapping.getNameScope()), insideFluxCorrection);
if (reactionParticipants[j] instanceof Product) {
if (rf.getFluxExpression().isZero()) {
rf.setFluxExpression(correctedReactionRate.flatten());
} else {
rf.setFluxExpression(Expression.add(rf.getFluxExpression(), correctedReactionRate.flatten()));
}
} else if (reactionParticipants[j] instanceof Reactant) {
if (rf.getFluxExpression().isZero()) {
rf.setFluxExpression(Expression.negate(correctedReactionRate).flatten());
} else {
rf.setFluxExpression(Expression.add(rf.getFluxExpression(), Expression.negate(correctedReactionRate).flatten()));
}
} else {
throw new RuntimeException("expected either FluxReactant or FluxProduct");
}
} else if (fr.getKinetics() instanceof LumpedKinetics) {
throw new RuntimeException("Lumped Kinetics for fluxes not yet supported");
} else {
throw new RuntimeException("unexpected Kinetic type in MembraneStructureAnalyzer.refreshResolvedFluxes()");
}
rf.getFluxExpression().bindExpression(mathMapping);
}
}
}
}
//
for (int i = 0; i < reactionSpecs.length; i++) {
if (reactionSpecs[i].isExcluded()) {
continue;
}
ReactionStep rs = reactionSpecs[i].getReactionStep();
if (rs.getStructure() != null && geoContext.getStructureMapping(rs.getStructure()).getGeometryClass() == surfaceClass) {
if (rs instanceof SimpleReaction) {
SimpleReaction sr = (SimpleReaction) rs;
ReactionParticipant[] rp_Array = sr.getReactionParticipants();
for (int k = 0; k < rp_Array.length; k++) {
if (rp_Array[k] instanceof Reactant || rp_Array[k] instanceof Product) {
SpeciesContextSpec rpSCS = mathMapping.getSimulationContext().getReactionContext().getSpeciesContextSpec(rp_Array[k].getSpeciesContext());
StructureMapping rpSM = mathMapping.getSimulationContext().getGeometryContext().getStructureMapping(rp_Array[k].getStructure());
if (rpSM.getGeometryClass() instanceof SubVolume && !rpSCS.isConstant()) {
//
// get ResolvedFlux for this species
//
ResolvedFlux rf = null;
for (int j = 0; j < resolvedFluxList.size(); j++) {
ResolvedFlux rf_tmp = (ResolvedFlux) resolvedFluxList.elementAt(j);
if (rf_tmp.getSpeciesContext() == rp_Array[k].getSpeciesContext()) {
rf = rf_tmp;
}
}
if (rf == null) {
VCUnitDefinition speciesFluxUnit = rp_Array[k].getSpeciesContext().getUnitDefinition().multiplyBy(unitSystem.getLengthUnit()).divideBy(unitSystem.getTimeUnit());
rf = new ResolvedFlux(rp_Array[k].getSpeciesContext(), speciesFluxUnit);
resolvedFluxList.addElement(rf);
}
if (rpSM.getGeometryClass() instanceof SubVolume && surfaceClass.isAdjacentTo((SubVolume) rpSM.getGeometryClass())) {
//
// for binding on inside or outside, add to ResolvedFlux.flux
//
Expression fluxRateExpression = getCorrectedRateExpression(sr, rp_Array[k], RateType.ResolvedFluxRate).renameBoundSymbols(mathMapping.getNameScope());
if (rf.getFluxExpression().isZero()) {
rf.setFluxExpression(fluxRateExpression);
} else {
rf.setFluxExpression(Expression.add(rf.getFluxExpression(), fluxRateExpression));
}
rf.getFluxExpression().bindExpression(mathMapping);
} else {
String structureName = ((rs.getStructure() != null) ? (rs.getStructure().getName()) : ("<null>"));
throw new Exception("In Application '" + mathMapping.getSimulationContext().getName() + "', SpeciesContext '" + rp_Array[k].getSpeciesContext().getName() + "' is not mapped adjacent to structure '" + structureName + "' but reacts there");
}
}
}
}
}
}
}
//
if (resolvedFluxList.size() > 0) {
resolvedFluxes = new ResolvedFlux[resolvedFluxList.size()];
resolvedFluxList.copyInto(resolvedFluxes);
} else {
resolvedFluxes = null;
}
}
use of cbit.vcell.model.ModelUnitSystem in project vcell by virtualcell.
the class ParticleMathMapping method combineHybrid.
private void combineHybrid() throws MappingException, ExpressionException, MatrixException, MathException, ModelException {
ArrayList<SpeciesContext> continuousSpecies = new ArrayList<SpeciesContext>();
ArrayList<ParticleVariable> continuousSpeciesParticleVars = new ArrayList<ParticleVariable>();
ArrayList<SpeciesContext> stochSpecies = new ArrayList<SpeciesContext>();
//
// categorize speciesContexts as continuous and stochastic
//
SpeciesContextSpec[] scsArray = getSimulationContext().getReactionContext().getSpeciesContextSpecs();
continuousSpecies = new ArrayList<SpeciesContext>();
stochSpecies = new ArrayList<SpeciesContext>();
for (SpeciesContextSpec speciesContextSpec : scsArray) {
if (!getSimulationContext().isStoch() || speciesContextSpec.isForceContinuous()) {
continuousSpecies.add(speciesContextSpec.getSpeciesContext());
Variable variable = getMathSymbolMapping().getVariable(speciesContextSpec.getSpeciesContext());
if (variable instanceof ParticleVariable) {
continuousSpeciesParticleVars.add((ParticleVariable) variable);
}
} else {
stochSpecies.add(speciesContextSpec.getSpeciesContext());
}
}
if (continuousSpecies.isEmpty()) {
return;
}
//
// create continuous mathDescription ... add stochastic variables and processes to the continuous Math and use this.
//
DiffEquMathMapping mathMapping = new DiffEquMathMapping(getSimulationContext(), callback, networkGenerationRequirements);
mathMapping.refresh(null);
MathDescription contMathDesc = mathMapping.getMathDescription();
//
// get list of all continuous variables
//
HashMap<String, Variable> allContinuousVars = new HashMap<String, Variable>();
Enumeration<Variable> enumVar = contMathDesc.getVariables();
while (enumVar.hasMoreElements()) {
Variable var = enumVar.nextElement();
allContinuousVars.put(var.getName(), var);
}
//
// replace those continuous variables and equations for stochastic speciesContexts
// with the particleVariables and particleProperties
// (ParticleJumpProcesses removed later)
//
ModelUnitSystem unitSystem = getSimulationContext().getModel().getUnitSystem();
for (SpeciesContext stochSpeciesContext : stochSpecies) {
Variable contVar = mathMapping.getMathSymbolMapping().getVariable(stochSpeciesContext);
Variable stochVar = getMathSymbolMapping().getVariable(stochSpeciesContext);
allContinuousVars.put(stochVar.getName(), stochVar);
//
// replace continuous "concentration" VolVariable/MemVariable for this particle with a Function for concentration
//
allContinuousVars.remove(contVar);
VCUnitDefinition sizeUnit = unitSystem.getLengthUnit().raiseTo(new RationalNumber(stochSpeciesContext.getStructure().getDimension()));
VCUnitDefinition stochasticDensityUnit = unitSystem.getStochasticSubstanceUnit().divideBy(sizeUnit);
VCUnitDefinition continuousDensityUnit = unitSystem.getConcentrationUnit(stochSpeciesContext.getStructure());
if (stochasticDensityUnit.isEquivalent(continuousDensityUnit)) {
allContinuousVars.put(contVar.getName(), new Function(contVar.getName(), new Expression(stochVar, getNameScope()), contVar.getDomain()));
} else {
Expression conversionFactorExp = getUnitFactor(continuousDensityUnit.divideBy(stochasticDensityUnit));
allContinuousVars.put(contVar.getName(), new Function(contVar.getName(), Expression.mult(new Expression(stochVar, getNameScope()), conversionFactorExp), contVar.getDomain()));
}
//
// remove continuous equation
//
Enumeration<SubDomain> contSubDomains = contMathDesc.getSubDomains();
while (contSubDomains.hasMoreElements()) {
SubDomain contSubDomain = contSubDomains.nextElement();
contSubDomain.removeEquation(contVar);
if (contSubDomain instanceof MembraneSubDomain) {
((MembraneSubDomain) contSubDomain).removeJumpCondition(contVar);
}
}
//
// remove all continuous variables for speciesContextSpec parameters (e.g. initial conditions, diffusion rates, boundary conditions, velocities)
//
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(stochSpeciesContext);
Parameter[] scsParameters = scs.getParameters();
for (Parameter parameter : scsParameters) {
Variable continuousScsParmVariable = mathMapping.getMathSymbolMapping().getVariable(parameter);
allContinuousVars.remove(continuousScsParmVariable);
}
//
// copy ParticleJumpProcess and ParticleProperties to the continuous math
//
SubDomain contSubDomain = contMathDesc.getSubDomain(contVar.getDomain().getName());
SubDomain stochSubDomain = mathDesc.getSubDomain(stochVar.getDomain().getName());
ParticleProperties particleProperties = stochSubDomain.getParticleProperties(stochVar);
contSubDomain.addParticleProperties(particleProperties);
}
//
// add all ParticleJumpProcesses to the continuous model
//
Enumeration<SubDomain> enumStochSubdomains = mathDesc.getSubDomains();
while (enumStochSubdomains.hasMoreElements()) {
SubDomain stochSubdomain = enumStochSubdomains.nextElement();
SubDomain contSubdomain = contMathDesc.getSubDomain(stochSubdomain.getName());
for (ParticleJumpProcess particleJumpProcess : stochSubdomain.getParticleJumpProcesses()) {
//
// modify "selection list" (particleVariables), probability rate, and actions if referenced particleVariable is to be "forced continuous"
//
ParticleVariable[] selectedParticles = particleJumpProcess.getParticleVariables();
for (ParticleVariable particleVariable : selectedParticles) {
if (continuousSpeciesParticleVars.contains(particleVariable)) {
particleJumpProcess.remove(particleVariable);
JumpProcessRateDefinition jumpProcessRateDefinition = particleJumpProcess.getParticleRateDefinition();
if (jumpProcessRateDefinition instanceof MacroscopicRateConstant) {
MacroscopicRateConstant macroscopicRateConstant = (MacroscopicRateConstant) jumpProcessRateDefinition;
macroscopicRateConstant.setExpression(Expression.mult(macroscopicRateConstant.getExpression(), new Expression(particleVariable, null)));
} else if (jumpProcessRateDefinition instanceof InteractionRadius) {
throw new MappingException("cannot adjust interaction radius for reaction process " + particleJumpProcess.getName() + ", particle " + particleVariable.getName() + " is continuous");
} else {
throw new MappingException("rate definition type " + jumpProcessRateDefinition.getClass().getSimpleName() + " not yet implemented for hybrid PDE/Particle math generation");
}
}
Iterator<Action> iterAction = particleJumpProcess.getActions().iterator();
while (iterAction.hasNext()) {
Action action = iterAction.next();
if (continuousSpeciesParticleVars.contains(action.getVar())) {
iterAction.remove();
}
}
}
if (!particleJumpProcess.getActions().isEmpty()) {
contSubdomain.addParticleJumpProcess(particleJumpProcess);
}
}
}
//
for (MathMappingParameter mathMappingParameter : fieldMathMappingParameters) {
if (mathMappingParameter instanceof UnitFactorParameter) {
String name = mathMappingParameter.getName();
if (!allContinuousVars.containsKey(name)) {
allContinuousVars.put(name, newFunctionOrConstant(name, mathMappingParameter.getExpression(), null));
}
}
}
//
// add constants and functions from the particle math that aren't already defined in the continuous math
//
Enumeration<Variable> enumVars = mathDesc.getVariables();
while (enumVars.hasMoreElements()) {
Variable var = enumVars.nextElement();
if (var instanceof Constant || var instanceof Function) {
String name = var.getName();
if (!allContinuousVars.containsKey(name)) {
allContinuousVars.put(name, var);
}
}
}
contMathDesc.setAllVariables(allContinuousVars.values().toArray(new Variable[0]));
mathDesc = contMathDesc;
//
for (int i = 0; i < fieldMathMappingParameters.length; i++) {
if (fieldMathMappingParameters[i] instanceof UnitFactorParameter) {
GeometryClass geometryClass = fieldMathMappingParameters[i].getGeometryClass();
Variable variable = newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass());
if (mathDesc.getVariable(variable.getName()) == null) {
mathDesc.addVariable(variable);
}
}
}
if (!mathDesc.isValid()) {
System.out.println(mathDesc.getVCML_database());
throw new MappingException("generated an invalid mathDescription: " + mathDesc.getWarning());
}
System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string begin ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
System.out.println(mathDesc.getVCML());
System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string end ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
}
use of cbit.vcell.model.ModelUnitSystem in project vcell by virtualcell.
the class RulebasedMathMapping method addStrictMassActionParticleJumpProcess.
private void addStrictMassActionParticleJumpProcess(VariableHash varHash, GeometryClass geometryClass, SubDomain subDomain, ReactionRule reactionRule, String jpName, ArrayList<ParticleVariable> reactantParticles, ArrayList<ParticleVariable> productParticles, ArrayList<Action> forwardActions, ArrayList<Action> reverseActions) throws ExpressionException, ExpressionBindingException, PropertyVetoException, MathException, MappingException {
String reactionRuleName = reactionRule.getName();
RbmKineticLaw kinetics = reactionRule.getKineticLaw();
RulebasedTransformation ruleBasedTransformation = ((RulebasedTransformation) getTransformation());
if (kinetics.getRateLawType() != RbmKineticLaw.RateLawType.MassAction) {
throw new RuntimeException("expecting mass action kinetics for reaction rule " + reactionRuleName);
}
//
// construct stochastic forward or reverse rate expression (separately). Transform from
// original expression of "concentrationRate" in terms of rateParameter and reactants/products in concentrations
// to
// new stochastic expression of "molecularRate" in terms of forwardRateParameter, reactants/products in molecules, structure sizes, and unit conversions.
//
// (1) concentrationRate = K * [s0] * [s1] [uM.s-1] or [molecules.um-3.s-1] or [molecules.um-2.s-1] (or other)
// (2) molecularRate = P * <s0> * <s1> [molecules.s-1]
//
// in this math description, we are using <s_i> [molecules], but original kinetics were in [s_i] [uM or molecules.um-2].
// so through a change in variable to get things in terms of <s_i>. <<<< Here P is the desired stochastic rate coefficient. >>>
//
// (3) let [s_i] = <s_i>/structsize(s_i)*unitConversionFactor(substanceunit([s_i])/substanceunit(<s_i>))
//
// in addition to the change in variables, we need to transform the entire expression from concentration/time to molecules/time
//
// (4) let molecularRate = concentrationRate * structSize(reaction) * unitConversionFactor(substanceunit(molecularRate)/substanceunit(concentrationRate))
//
// (5) in general, concentationRate = K * PRODUCT([s_i])
//
// change of variables into stochastic variables used in MathDescription, substituting (3) into (5)
//
// (6) concentrationRate = K * PRODUCT(<s_i>/structsize(s_i)*unitConversionFactor(substanceunit([s_i])/substanceunit(<s_i>)))
//
// reordering to separate the sizes, the unit conversions and the <s_i>
//
// (7) concentrationRate = K * PRODUCT(<s_i>) * PRODUCT(1/structsize(s_i)) * unitConversionFactor(PRODUCT(substanceunit([s_i])/substanceunit(<s_i>)))
//
// combining (4) and (7)
//
// (8) molecularRate = K * PRODUCT(<s_i>) * PRODUCT(1/structsize(s_i)) * unitConversionFactor(PRODUCT(substanceunit([s_i])/substanceunit(<s_i>))) * structSize(reaction) * unitConversionFactor(substanceunit(molecularRate)/substanceunit(concentrationRate))
//
// collecting terms of sizes and unit conversions
//
// (9) molecularRate = K * PRODUCT(<s_i>) * structSize(reaction) / PRODUCT(structsize(s_i)) * unitConversionFactor(substanceunit(molecularRate)/substanceunit(concentrationRate) * PRODUCT(substanceunit([s_i])/substanceunit(<s_i>)))
//
// (10) molecularRate = K * PRODUCT(<s_i>) * sizeFactor * unitConversionFactor(substanceConversionUnit)
//
// where
//
// (11) sizeFactor = structSize(reaction) / PRODUCT(structsize(s_i))
// (12) substanceConversionUnit = substanceunit(molecularRate)/substanceunit(concentrationRate) * PRODUCT(substanceunit([s_i])/substanceunit(<s_i>))
//
// The ParticleJumpCondition wants a single new rate stochastic, P from equation (2). Note that PRODUCT(<s_i>) will be captured separately the the reactantPatterns.
// comparing (2) and (10) we have found P.
//
// (13) P = K * sizeFactor * unitConversionFactor(substanceConversionUnit)
//
// the framework also needs the proper unit for P
//
// (14) Unit(P) = Unit(K) * Unit(sizeFactor) * substanceConversionUnit
//
//
ModelUnitSystem modelUnitSystem = getSimulationContext().getModel().getUnitSystem();
VCUnitDefinition stochasticSubstanceUnit = modelUnitSystem.getStochasticSubstanceUnit();
VCUnitDefinition reactionRuleSubstanceUnit = modelUnitSystem.getSubstanceUnit(reactionRule.getStructure());
int forwardRuleIndex = 0;
//
// get forward rate parameter and make sure it is constant valued.
//
Parameter forward_rateParameter = kinetics.getLocalParameter(RbmKineticLawParameterType.MassActionForwardRate);
Expression substitutedForwardRate = MathUtilities.substituteModelParameters(forward_rateParameter.getExpression(), reactionRule.getNameScope().getScopedSymbolTable());
if (!substitutedForwardRate.flatten().isNumeric()) {
throw new MappingException("forward rate constant for reaction rule " + reactionRule.getName() + " is not constant");
}
//
// create forward sizeExp and forward unitFactor
//
VCUnitDefinition forward_substanceConversionUnit = stochasticSubstanceUnit.divideBy(reactionRuleSubstanceUnit);
VCUnitDefinition forward_sizeFactorUnit = reactionRule.getStructure().getStructureSize().getUnitDefinition();
Expression forward_sizeFactor = new Expression(reactionRule.getStructure().getStructureSize(), getNameScope());
for (ReactantPattern reactantPattern : reactionRule.getReactantPatterns()) {
Expression reactantSizeExp = new Expression(reactantPattern.getStructure().getStructureSize(), getNameScope());
VCUnitDefinition reactantSizeUnit = reactantPattern.getStructure().getStructureSize().getUnitDefinition();
VCUnitDefinition reactantSubstanceUnit = modelUnitSystem.getSubstanceUnit(reactantPattern.getStructure());
forward_sizeFactor = Expression.div(forward_sizeFactor, reactantSizeExp);
forward_sizeFactorUnit = forward_sizeFactorUnit.divideBy(reactantSizeUnit);
forward_substanceConversionUnit = forward_substanceConversionUnit.multiplyBy(reactantSubstanceUnit).divideBy(stochasticSubstanceUnit);
}
// simplify sizeFactor (often has size/size/size)
try {
forward_sizeFactor = RationalExpUtils.getRationalExp(forward_sizeFactor).simplifyAsExpression();
forward_sizeFactor.bindExpression(getSimulationContext().getModel());
} catch (ParseException e) {
e.printStackTrace();
}
Expression forward_rateExp = Expression.mult(new Expression(forward_rateParameter, getNameScope()), forward_sizeFactor, getUnitFactor(forward_substanceConversionUnit)).flatten();
VCUnitDefinition forward_rateUnit = forward_rateParameter.getUnitDefinition().multiplyBy(forward_sizeFactorUnit).multiplyBy(forward_substanceConversionUnit);
ProbabilityParameter forward_probParm = addProbabilityParameter(PARAMETER_PROBABILITYRATE_PREFIX + jpName, forward_rateExp, PARAMETER_ROLE_P, forward_rateUnit, reactionRule);
// add probability to function or constant
varHash.addVariable(newFunctionOrConstant(getMathSymbol(forward_probParm, geometryClass), getIdentifierSubstitutions(forward_rateExp, forward_rateUnit, geometryClass), geometryClass));
// add forward ParticleJumpProcess
String forward_name = reactionRuleName;
Expression forward_rate = getIdentifierSubstitutions(new Expression(forward_probParm, getNameScope()), forward_probParm.getUnitDefinition(), geometryClass);
JumpProcessRateDefinition forward_rateDefinition = new MacroscopicRateConstant(forward_rate);
ReactionRuleAnalysisReport rrarBiomodelForward = ruleBasedTransformation.getRulesForwardMap().get(reactionRule);
ProcessSymmetryFactor forwardSymmetryFactor = new ProcessSymmetryFactor(rrarBiomodelForward.getSymmetryFactor());
ParticleJumpProcess forward_particleJumpProcess = new ParticleJumpProcess(forward_name, reactantParticles, forward_rateDefinition, forwardActions, forwardSymmetryFactor);
subDomain.addParticleJumpProcess(forward_particleJumpProcess);
//
for (ReactionRule rr : getSimulationContext().getModel().getRbmModelContainer().getReactionRuleList()) {
if (rr == reactionRule) {
break;
}
forwardRuleIndex++;
if (rr.isReversible()) {
forwardRuleIndex++;
}
}
//
if (reactionRule.isReversible()) {
Parameter reverse_rateParameter = kinetics.getLocalParameter(RbmKineticLawParameterType.MassActionReverseRate);
if (reverse_rateParameter == null || reverse_rateParameter.getExpression() == null) {
throw new MappingException("reverse rate constant for reaction rule " + reactionRule.getName() + " is missing");
}
{
Expression substitutedReverseRate = MathUtilities.substituteModelParameters(reverse_rateParameter.getExpression(), reactionRule.getNameScope().getScopedSymbolTable());
if (!substitutedReverseRate.flatten().isNumeric()) {
throw new MappingException("reverse rate constant for reaction rule " + reactionRule.getName() + " is not constant");
}
}
//
// create reverse sizeExp and reverse unitFactor
//
VCUnitDefinition reverse_substanceConversionUnit = stochasticSubstanceUnit.divideBy(reactionRuleSubstanceUnit);
VCUnitDefinition reverse_sizeFactorUnit = reactionRule.getStructure().getStructureSize().getUnitDefinition();
Expression reverse_sizeFactor = new Expression(reactionRule.getStructure().getStructureSize(), getNameScope());
for (ProductPattern productPattern : reactionRule.getProductPatterns()) {
Expression reactantSizeExp = new Expression(productPattern.getStructure().getStructureSize(), getNameScope());
VCUnitDefinition reactantSizeUnit = productPattern.getStructure().getStructureSize().getUnitDefinition();
VCUnitDefinition reactantSubstanceUnit = modelUnitSystem.getSubstanceUnit(productPattern.getStructure());
reverse_sizeFactor = Expression.div(reverse_sizeFactor, reactantSizeExp);
reverse_sizeFactorUnit = reverse_sizeFactorUnit.divideBy(reactantSizeUnit);
reverse_substanceConversionUnit = reverse_substanceConversionUnit.multiplyBy(reactantSubstanceUnit).divideBy(stochasticSubstanceUnit);
}
// simplify sizeFactor (often has size/size/size)
try {
reverse_sizeFactor = RationalExpUtils.getRationalExp(reverse_sizeFactor).simplifyAsExpression();
reverse_sizeFactor.bindExpression(getSimulationContext().getModel());
} catch (ParseException e) {
e.printStackTrace();
}
Expression reverse_rateExp = Expression.mult(new Expression(reverse_rateParameter, getNameScope()), reverse_sizeFactor, getUnitFactor(reverse_substanceConversionUnit)).flatten();
VCUnitDefinition reverse_rateUnit = reverse_rateParameter.getUnitDefinition().multiplyBy(reverse_sizeFactorUnit).multiplyBy(reverse_substanceConversionUnit);
// if the reaction has forward rate (Mass action,HMMs), or don't have either forward or reverse rate (some other rate laws--like general)
// we process it as forward reaction
// get jump process name
ProbabilityParameter reverse_probParm = addProbabilityParameter(PARAMETER_PROBABILITYRATE_PREFIX + jpName + "_reverse", reverse_rateExp, PARAMETER_ROLE_P_reverse, reverse_rateUnit, reactionRule);
// add probability to function or constant
varHash.addVariable(newFunctionOrConstant(getMathSymbol(reverse_probParm, geometryClass), getIdentifierSubstitutions(reverse_rateExp, reverse_rateUnit, geometryClass), geometryClass));
// add reverse ParticleJumpProcess
Expression reverse_rate = getIdentifierSubstitutions(new Expression(reverse_probParm, getNameScope()), reverse_probParm.getUnitDefinition(), geometryClass);
String reverse_name = reactionRuleName + "_reverse";
JumpProcessRateDefinition reverse_rateDefinition = new MacroscopicRateConstant(reverse_rate);
ReactionRuleAnalysisReport rrarBiomodelReverse = ruleBasedTransformation.getRulesReverseMap().get(reactionRule);
ProcessSymmetryFactor reverseSymmetryFactor = new ProcessSymmetryFactor(rrarBiomodelReverse.getSymmetryFactor());
ParticleJumpProcess reverse_particleJumpProcess = new ParticleJumpProcess(reverse_name, productParticles, reverse_rateDefinition, reverseActions, reverseSymmetryFactor);
subDomain.addParticleJumpProcess(reverse_particleJumpProcess);
//
// check reverse direction mapping and operations with RuleAnalysis.
//
int reverseRuleIndex = forwardRuleIndex + 1;
ReactionRuleAnalysisReport rrar = ruleBasedTransformation.getRulesReverseMap().get(reactionRule);
jumpProcessMap.put(reverse_particleJumpProcess, rrar);
}
}
Aggregations