Example 6 with Reactant
use of cbit.vcell.model.Reactant in project vcell by virtualcell.
the class ModelProcessEquation method parseReaction.
public static ReactionParticipant[] parseReaction(ReactionStep reactionStep, Model model, String equationString) throws ExpressionException, PropertyVetoException {
int gotoIndex = equationString.indexOf(REACTION_GOESTO);
if (gotoIndex < 1 && equationString.length() == 0) {
throw new ExpressionException("Syntax error! " + REACTION_GOESTO + " not found. (e.g. a+b->c)");
}
if (reactionStep == null) {
return null;
}
String leftHand = equationString.substring(0, gotoIndex);
String rightHand = equationString.substring(gotoIndex + REACTION_GOESTO.length());
StringTokenizer st = new StringTokenizer(leftHand, "+");
ArrayList<ReactionParticipant> rplist = new ArrayList<ReactionParticipant>();
HashMap<String, SpeciesContext> speciesContextMap = new HashMap<String, SpeciesContext>();
Structure rxnStructure = reactionStep.getStructure();
while (st.hasMoreElements()) {
String nextToken = st.nextToken().trim();
if (nextToken.length() == 0) {
continue;
}
int stoichiIndex = 0;
while (true) {
if (Character.isDigit(nextToken.charAt(stoichiIndex))) {
stoichiIndex++;
} else {
break;
}
}
int stoichi = 1;
String tmp = nextToken.substring(0, stoichiIndex);
if (tmp.length() > 0) {
stoichi = Integer.parseInt(tmp);
}
String var = nextToken.substring(stoichiIndex).trim();
SpeciesContext sc = model.getSpeciesContext(var);
if (sc == null) {
sc = speciesContextMap.get(var);
if (sc == null) {
Species species = model.getSpecies(var);
if (species == null) {
species = new Species(var, null);
}
sc = new SpeciesContext(species, rxnStructure);
sc.setName(var);
speciesContextMap.put(var, sc);
}
}
// if (reactionStep instanceof SimpleReaction) {
rplist.add(new Reactant(null, (SimpleReaction) reactionStep, sc, stoichi));
// } else if (reactionStep instanceof FluxReaction) {
// rplist.add(new Flux(null, (FluxReaction) reactionStep, sc));
// }
}
st = new StringTokenizer(rightHand, "+");
while (st.hasMoreElements()) {
String nextToken = st.nextToken().trim();
if (nextToken.length() == 0) {
continue;
}
int stoichiIndex = 0;
while (true) {
if (Character.isDigit(nextToken.charAt(stoichiIndex))) {
stoichiIndex++;
} else {
break;
}
}
int stoichi = 1;
String tmp = nextToken.substring(0, stoichiIndex);
if (tmp.length() > 0) {
stoichi = Integer.parseInt(tmp);
}
String var = nextToken.substring(stoichiIndex);
SpeciesContext sc = model.getSpeciesContext(var);
if (sc == null) {
sc = speciesContextMap.get(var);
if (sc == null) {
Species species = model.getSpecies(var);
if (species == null) {
species = new Species(var, null);
}
sc = new SpeciesContext(species, rxnStructure);
sc.setName(var);
speciesContextMap.put(var, sc);
}
}
// if (reactionStep instanceof SimpleReaction) {
rplist.add(new Product(null, (SimpleReaction) reactionStep, sc, stoichi));
// } else if (reactionStep instanceof FluxReaction) {
// rplist.add(new Flux(null, (FluxReaction) reactionStep, sc));
// }
}
return rplist.toArray(new ReactionParticipant[0]);
}
Also used :
SimpleReaction(cbit.vcell.model.SimpleReaction)
HashMap(java.util.HashMap)
ArrayList(java.util.ArrayList)
Product(cbit.vcell.model.Product)
SpeciesContext(cbit.vcell.model.SpeciesContext)
Reactant(cbit.vcell.model.Reactant)
ExpressionException(cbit.vcell.parser.ExpressionException)
StringTokenizer(java.util.StringTokenizer)
Structure(cbit.vcell.model.Structure)
ReactionParticipant(cbit.vcell.model.ReactionParticipant)
Species(cbit.vcell.model.Species)
Example 7 with Reactant
use of cbit.vcell.model.Reactant in project vcell by virtualcell.
the class ParticleMathMapping method refreshMathDescription.
/**
* This method was created in VisualAge.
*/
private void refreshMathDescription() throws MappingException, MatrixException, MathException, ExpressionException, ModelException {
getSimulationContext().checkValidity();
if (getSimulationContext().getGeometry().getDimension() == 0) {
throw new MappingException("particle math mapping requires spatial geometry - dimension >= 1");
}
StructureMapping[] structureMappings = getSimulationContext().getGeometryContext().getStructureMappings();
for (int i = 0; i < structureMappings.length; i++) {
if (structureMappings[i] instanceof MembraneMapping) {
if (((MembraneMapping) structureMappings[i]).getCalculateVoltage()) {
throw new MappingException("electric potential not yet supported for particle models");
}
}
}
//
// fail if any events
//
BioEvent[] bioEvents = getSimulationContext().getBioEvents();
if (bioEvents != null && bioEvents.length > 0) {
throw new MappingException("events not yet supported for particle-based models");
}
//
// gather only those reactionSteps that are not "excluded"
//
ReactionSpec[] reactionSpecs = getSimulationContext().getReactionContext().getReactionSpecs();
Vector<ReactionStep> rsList = new Vector<ReactionStep>();
for (int i = 0; i < reactionSpecs.length; i++) {
if (reactionSpecs[i].isExcluded() == false) {
if (reactionSpecs[i].isFast()) {
throw new MappingException("fast reactions not supported for particle models");
}
rsList.add(reactionSpecs[i].getReactionStep());
}
}
ReactionStep[] reactionSteps = new ReactionStep[rsList.size()];
rsList.copyInto(reactionSteps);
//
for (int i = 0; i < reactionSteps.length; i++) {
Kinetics.UnresolvedParameter[] unresolvedParameters = reactionSteps[i].getKinetics().getUnresolvedParameters();
if (unresolvedParameters != null && unresolvedParameters.length > 0) {
StringBuffer buffer = new StringBuffer();
for (int j = 0; j < unresolvedParameters.length; j++) {
if (j > 0) {
buffer.append(", ");
}
buffer.append(unresolvedParameters[j].getName());
}
throw new MappingException(reactionSteps[i].getDisplayType() + " '" + reactionSteps[i].getName() + "' contains unresolved identifier(s): " + buffer);
}
}
//
// temporarily place all variables in a hashtable (before binding) and discarding duplicates (check for equality)
//
VariableHash varHash = new VariableHash();
// //
// // verify that all structures are mapped to geometry classes and all geometry classes are mapped to a structure
// //
// Structure structures[] = getSimulationContext().getGeometryContext().getModel().getStructures();
// for (int i = 0; i < structures.length; i++){
// StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(structures[i]);
// if (sm==null || (sm.getGeometryClass() == null)){
// throw new MappingException("model structure '"+structures[i].getName()+"' not mapped to a geometry subdomain");
// }
// if (sm.getUnitSizeParameter()!=null){
// Expression unitSizeExp = sm.getUnitSizeParameter().getExpression();
// if(unitSizeExp != null)
// {
// try {
// double unitSize = unitSizeExp.evaluateConstant();
// if (unitSize != 1.0){
// throw new MappingException("model structure '"+sm.getStructure().getName()+"' unit size = "+unitSize+" != 1.0 ... partial volume or surface mapping not yet supported for particles");
// }
// }catch (ExpressionException e){
// e.printStackTrace(System.out);
// throw new MappingException("couldn't evaluate unit size for model structure '"+sm.getStructure().getName()+"' : "+e.getMessage());
// }
// }
// }
// }
// {
// GeometryClass[] geometryClass = getSimulationContext().getGeometryContext().getGeometry().getGeometryClasses();
// for (int i = 0; i < geometryClass.length; i++){
// Structure[] mappedStructures = getSimulationContext().getGeometryContext().getStructuresFromGeometryClass(geometryClass[i]);
// if (mappedStructures==null || mappedStructures.length==0){
// throw new MappingException("geometryClass '"+geometryClass[i].getName()+"' not mapped from a model structure");
// }
// }
// }
// deals with model parameters
Model model = getSimulationContext().getModel();
ModelUnitSystem modelUnitSystem = model.getUnitSystem();
ModelParameter[] modelParameters = model.getModelParameters();
// populate in globalParameterVariants hashtable
for (int j = 0; j < modelParameters.length; j++) {
Expression modelParamExpr = modelParameters[j].getExpression();
GeometryClass geometryClass = getDefaultGeometryClass(modelParamExpr);
modelParamExpr = getIdentifierSubstitutions(modelParamExpr, modelParameters[j].getUnitDefinition(), geometryClass);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(modelParameters[j], geometryClass), modelParamExpr, geometryClass));
}
//
// create new MathDescription (based on simContext's previous MathDescription if possible)
//
MathDescription oldMathDesc = getSimulationContext().getMathDescription();
mathDesc = null;
if (oldMathDesc != null) {
if (oldMathDesc.getVersion() != null) {
mathDesc = new MathDescription(oldMathDesc.getVersion());
} else {
mathDesc = new MathDescription(oldMathDesc.getName());
}
} else {
mathDesc = new MathDescription(getSimulationContext().getName() + "_generated");
}
//
// volume particle variables
//
Enumeration<SpeciesContextMapping> enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = enum1.nextElement();
if (scm.getVariable() instanceof ParticleVariable) {
if (!(mathDesc.getVariable(scm.getVariable().getName()) instanceof ParticleVariable)) {
varHash.addVariable(scm.getVariable());
}
}
}
varHash.addVariable(new Constant(getMathSymbol(model.getPI_CONSTANT(), null), getIdentifierSubstitutions(model.getPI_CONSTANT().getExpression(), model.getPI_CONSTANT().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getFARADAY_CONSTANT(), null), getIdentifierSubstitutions(model.getFARADAY_CONSTANT().getExpression(), model.getFARADAY_CONSTANT().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getFARADAY_CONSTANT_NMOLE(), null), getIdentifierSubstitutions(model.getFARADAY_CONSTANT_NMOLE().getExpression(), model.getFARADAY_CONSTANT_NMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getGAS_CONSTANT(), null), getIdentifierSubstitutions(model.getGAS_CONSTANT().getExpression(), model.getGAS_CONSTANT().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getTEMPERATURE(), null), getIdentifierSubstitutions(new Expression(getSimulationContext().getTemperatureKelvin()), model.getTEMPERATURE().getUnitDefinition(), null)));
//
for (int j = 0; j < structureMappings.length; j++) {
if (structureMappings[j] instanceof MembraneMapping) {
MembraneMapping membraneMapping = (MembraneMapping) structureMappings[j];
GeometryClass geometryClass = membraneMapping.getGeometryClass();
//
// don't calculate voltage, still may need it though
//
Parameter initialVoltageParm = membraneMapping.getInitialVoltageParameter();
Variable voltageFunction = newFunctionOrConstant(getMathSymbol(membraneMapping.getMembrane().getMembraneVoltage(), geometryClass), getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), geometryClass), geometryClass);
varHash.addVariable(voltageFunction);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(membraneMapping.getMembrane().getMembraneVoltage(), membraneMapping.getGeometryClass()), getIdentifierSubstitutions(membraneMapping.getInitialVoltageParameter().getExpression(), membraneMapping.getInitialVoltageParameter().getUnitDefinition(), membraneMapping.getGeometryClass()), membraneMapping.getGeometryClass()));
}
}
//
for (int j = 0; j < reactionSteps.length; j++) {
ReactionStep rs = reactionSteps[j];
if (getSimulationContext().getReactionContext().getReactionSpec(rs).isExcluded()) {
continue;
}
Kinetics.KineticsParameter[] parameters = rs.getKinetics().getKineticsParameters();
GeometryClass geometryClass = null;
if (rs.getStructure() != null) {
geometryClass = getSimulationContext().getGeometryContext().getStructureMapping(rs.getStructure()).getGeometryClass();
}
if (parameters != null) {
for (int i = 0; i < parameters.length; i++) {
// Reaction rate, currentDensity, LumpedCurrent and null parameters are not going to displayed in the particle math description.
if (((parameters[i].getRole() == Kinetics.ROLE_CurrentDensity) || (parameters[i].getRole() == Kinetics.ROLE_LumpedCurrent) || (parameters[i].getRole() == Kinetics.ROLE_ReactionRate)) || (parameters[i].getExpression() == null)) {
continue;
}
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameters[i], geometryClass), getIdentifierSubstitutions(parameters[i].getExpression(), parameters[i].getUnitDefinition(), geometryClass), geometryClass));
}
}
}
//
// initial constants (either function or constant)
//
SpeciesContextSpec[] speciesContextSpecs = getSimulationContext().getReactionContext().getSpeciesContextSpecs();
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpecParameter initParm = null;
Expression initExpr = null;
if (getSimulationContext().isUsingConcentration()) {
initParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
initExpr = new Expression(initParm.getExpression());
// if (speciesContextSpecs[i].getSpeciesContext().getStructure() instanceof Feature) {
// initExpr = Expression.div(initExpr, new Expression(model.getKMOLE, getNameScope())).flatten();
// }
} else {
initParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_InitialCount);
initExpr = new Expression(initParm.getExpression());
}
if (initExpr != null) {
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
String[] symbols = initExpr.getSymbols();
// Check if 'initExpr' has other speciesContexts in its expression, need to replace it with 'spContext_init'
for (int j = 0; symbols != null && j < symbols.length; j++) {
// if symbol is a speciesContext, replacing it with a reference to initial condition for that speciesContext.
SpeciesContext spC = null;
SymbolTableEntry ste = initExpr.getSymbolBinding(symbols[j]);
if (ste instanceof SpeciesContextSpecProxyParameter) {
SpeciesContextSpecProxyParameter spspp = (SpeciesContextSpecProxyParameter) ste;
if (spspp.getTarget() instanceof SpeciesContext) {
spC = (SpeciesContext) spspp.getTarget();
SpeciesContextSpec spcspec = getSimulationContext().getReactionContext().getSpeciesContextSpec(spC);
SpeciesContextSpecParameter spCInitParm = spcspec.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
// if initConc param expression is null, try initCount
if (spCInitParm.getExpression() == null) {
spCInitParm = spcspec.getParameterFromRole(SpeciesContextSpec.ROLE_InitialCount);
}
// need to get init condn expression, but can't get it from getMathSymbol() (mapping between bio and math), hence get it as below.
Expression scsInitExpr = new Expression(spCInitParm, getNameScope());
// scsInitExpr.bindExpression(this);
initExpr.substituteInPlace(new Expression(spC.getName()), scsInitExpr);
}
}
}
// now create the appropriate function for the current speciesContextSpec.
varHash.addVariable(newFunctionOrConstant(getMathSymbol(initParm, sm.getGeometryClass()), getIdentifierSubstitutions(initExpr, initParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
//
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpec.SpeciesContextSpecParameter diffParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_DiffusionRate);
if (diffParm != null) {
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(diffParm, sm.getGeometryClass()), getIdentifierSubstitutions(diffParm.getExpression(), diffParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
//
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpec.SpeciesContextSpecParameter bc_xm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueXm);
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
if (bc_xm != null && (bc_xm.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_xm, sm.getGeometryClass()), getIdentifierSubstitutions(bc_xm.getExpression(), bc_xm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_xp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueXp);
if (bc_xp != null && (bc_xp.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_xp, sm.getGeometryClass()), getIdentifierSubstitutions(bc_xp.getExpression(), bc_xp.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_ym = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueYm);
if (bc_ym != null && (bc_ym.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_ym, sm.getGeometryClass()), getIdentifierSubstitutions(bc_ym.getExpression(), bc_ym.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_yp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueYp);
if (bc_yp != null && (bc_yp.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_yp, sm.getGeometryClass()), getIdentifierSubstitutions(bc_yp.getExpression(), bc_yp.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_zm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueZm);
if (bc_zm != null && (bc_zm.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_zm, sm.getGeometryClass()), getIdentifierSubstitutions(bc_zm.getExpression(), bc_zm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_zp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueZp);
if (bc_zp != null && (bc_zp.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_zp, sm.getGeometryClass()), getIdentifierSubstitutions(bc_zp.getExpression(), bc_zp.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
//
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpec.SpeciesContextSpecParameter advection_velX = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityX);
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
GeometryClass geometryClass = sm.getGeometryClass();
if (advection_velX != null && (advection_velX.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velX, geometryClass), getIdentifierSubstitutions(advection_velX.getExpression(), advection_velX.getUnitDefinition(), geometryClass), geometryClass));
}
SpeciesContextSpec.SpeciesContextSpecParameter advection_velY = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityY);
if (advection_velY != null && (advection_velY.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velY, geometryClass), getIdentifierSubstitutions(advection_velY.getExpression(), advection_velY.getUnitDefinition(), geometryClass), geometryClass));
}
SpeciesContextSpec.SpeciesContextSpecParameter advection_velZ = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityZ);
if (advection_velZ != null && (advection_velZ.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velZ, geometryClass), getIdentifierSubstitutions(advection_velZ.getExpression(), advection_velZ.getUnitDefinition(), geometryClass), geometryClass));
}
}
//
// constant species (either function or constant)
//
enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
if (scm.getVariable() instanceof Constant) {
varHash.addVariable(scm.getVariable());
}
}
//
// conversion factors
//
varHash.addVariable(new Constant(getMathSymbol(model.getKMOLE(), null), getIdentifierSubstitutions(model.getKMOLE().getExpression(), model.getKMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getN_PMOLE(), null), getIdentifierSubstitutions(model.getN_PMOLE().getExpression(), model.getN_PMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getKMILLIVOLTS(), null), getIdentifierSubstitutions(model.getKMILLIVOLTS().getExpression(), model.getKMILLIVOLTS().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getK_GHK(), null), getIdentifierSubstitutions(model.getK_GHK().getExpression(), model.getK_GHK().getUnitDefinition(), null)));
//
for (int i = 0; i < structureMappings.length; i++) {
StructureMapping sm = structureMappings[i];
if (getSimulationContext().getGeometry().getDimension() == 0) {
StructureMappingParameter sizeParm = sm.getSizeParameter();
if (sizeParm != null && sizeParm.getExpression() != null) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(sizeParm, sm.getGeometryClass()), getIdentifierSubstitutions(sizeParm.getExpression(), sizeParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
} else {
if (sm instanceof MembraneMapping) {
MembraneMapping mm = (MembraneMapping) sm;
StructureMappingParameter volFrac = mm.getVolumeFractionParameter();
if (volFrac != null && volFrac.getExpression() != null) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(volFrac, sm.getGeometryClass()), getIdentifierSubstitutions(volFrac.getExpression(), volFrac.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
StructureMappingParameter surfToVol = mm.getSurfaceToVolumeParameter();
if (surfToVol != null && surfToVol.getExpression() != null) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(surfToVol, sm.getGeometryClass()), getIdentifierSubstitutions(surfToVol.getExpression(), surfToVol.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
}
} else {
Parameter parm = sm.getParameterFromRole(StructureMapping.ROLE_AreaPerUnitArea);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SurfaceClass) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_AreaPerUnitVolume);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_VolumePerUnitArea);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SurfaceClass) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_VolumePerUnitVolume);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_Size);
if (parm != null && parm.getExpression() != null) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
}
//
// functions
//
enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
if (scm.getVariable() == null && scm.getDependencyExpression() != null) {
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(scm.getSpeciesContext().getStructure());
Variable dependentVariable = newFunctionOrConstant(getMathSymbol(scm.getSpeciesContext(), sm.getGeometryClass()), getIdentifierSubstitutions(scm.getDependencyExpression(), scm.getSpeciesContext().getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass());
dependentVariable.setDomain(new Domain(sm.getGeometryClass()));
varHash.addVariable(dependentVariable);
}
}
//
for (int i = 0; i < fieldMathMappingParameters.length; i++) {
if (fieldMathMappingParameters[i] instanceof UnitFactorParameter) {
GeometryClass geometryClass = fieldMathMappingParameters[i].getGeometryClass();
varHash.addVariable(newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass()));
}
}
//
// set Variables to MathDescription all at once with the order resolved by "VariableHash"
//
mathDesc.setAllVariables(varHash.getAlphabeticallyOrderedVariables());
//
if (getSimulationContext().getGeometryContext().getGeometry() != null) {
try {
mathDesc.setGeometry(getSimulationContext().getGeometryContext().getGeometry());
} catch (java.beans.PropertyVetoException e) {
e.printStackTrace(System.out);
throw new MappingException("failure setting geometry " + e.getMessage());
}
} else {
throw new MappingException("geometry must be defined");
}
//
// create subdomains (volume and surfaces)
//
GeometryClass[] geometryClasses = getSimulationContext().getGeometryContext().getGeometry().getGeometryClasses();
for (int k = 0; k < geometryClasses.length; k++) {
if (geometryClasses[k] instanceof SubVolume) {
SubVolume subVolume = (SubVolume) geometryClasses[k];
//
// get priority of subDomain
//
// now does not have to match spatial feature, *BUT* needs to be unique
int priority = k;
//
// create subDomain
//
CompartmentSubDomain subDomain = new CompartmentSubDomain(subVolume.getName(), priority);
mathDesc.addSubDomain(subDomain);
//
// assign boundary condition types
//
StructureMapping[] mappedSMs = getSimulationContext().getGeometryContext().getStructureMappings(subVolume);
FeatureMapping mappedFM = null;
for (int i = 0; i < mappedSMs.length; i++) {
if (mappedSMs[i] instanceof FeatureMapping) {
if (mappedFM != null) {
lg.warn("WARNING:::: MathMapping.refreshMathDescription() ... assigning boundary condition types not unique");
}
mappedFM = (FeatureMapping) mappedSMs[i];
}
}
if (mappedFM != null) {
subDomain.setBoundaryConditionXm(mappedFM.getBoundaryConditionTypeXm());
subDomain.setBoundaryConditionXp(mappedFM.getBoundaryConditionTypeXp());
if (getSimulationContext().getGeometry().getDimension() > 1) {
subDomain.setBoundaryConditionYm(mappedFM.getBoundaryConditionTypeYm());
subDomain.setBoundaryConditionYp(mappedFM.getBoundaryConditionTypeYp());
}
if (getSimulationContext().getGeometry().getDimension() > 2) {
subDomain.setBoundaryConditionZm(mappedFM.getBoundaryConditionTypeZm());
subDomain.setBoundaryConditionZp(mappedFM.getBoundaryConditionTypeZp());
}
}
} else if (geometryClasses[k] instanceof SurfaceClass) {
SurfaceClass surfaceClass = (SurfaceClass) geometryClasses[k];
// determine membrane inside and outside subvolume
// this preserves backward compatibility so that membrane subdomain
// inside and outside correspond to structure hierarchy when present
Pair<SubVolume, SubVolume> ret = DiffEquMathMapping.computeBoundaryConditionSource(model, simContext, surfaceClass);
SubVolume innerSubVolume = ret.one;
SubVolume outerSubVolume = ret.two;
//
// create subDomain
//
CompartmentSubDomain outerCompartment = mathDesc.getCompartmentSubDomain(outerSubVolume.getName());
CompartmentSubDomain innerCompartment = mathDesc.getCompartmentSubDomain(innerSubVolume.getName());
MembraneSubDomain memSubDomain = new MembraneSubDomain(innerCompartment, outerCompartment, surfaceClass.getName());
mathDesc.addSubDomain(memSubDomain);
}
}
//
// create Particle Contexts for all Particle Variables
//
Enumeration<SpeciesContextMapping> enumSCM = getSpeciesContextMappings();
Expression unitFactor = getUnitFactor(modelUnitSystem.getStochasticSubstanceUnit().divideBy(modelUnitSystem.getVolumeSubstanceUnit()));
while (enumSCM.hasMoreElements()) {
SpeciesContextMapping scm = enumSCM.nextElement();
SpeciesContext sc = scm.getSpeciesContext();
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure());
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
if (scm.getVariable() instanceof ParticleVariable && scm.getDependencyExpression() == null) {
ParticleVariable particleVariable = (ParticleVariable) scm.getVariable();
//
// initial distribution of particles
//
ArrayList<ParticleInitialCondition> particleInitialConditions = new ArrayList<ParticleInitialCondition>();
ParticleInitialCondition pic = null;
if (getSimulationContext().isUsingConcentration()) {
Expression initialDistribution = scs.getInitialConcentrationParameter().getExpression() == null ? null : new Expression(getMathSymbol(scs.getInitialConcentrationParameter(), sm.getGeometryClass()));
if (particleVariable instanceof VolumeParticleVariable) {
initialDistribution = Expression.mult(initialDistribution, unitFactor);
}
pic = new ParticleInitialConditionConcentration(initialDistribution);
} else {
Expression initialCount = scs.getInitialCountParameter().getExpression() == null ? null : new Expression(getMathSymbol(scs.getInitialCountParameter(), sm.getGeometryClass()));
if (initialCount == null) {
throw new MappingException("initialCount not defined for speciesContext " + scs.getSpeciesContext().getName());
}
Expression locationX = new Expression("u");
Expression locationY = new Expression("u");
Expression locationZ = new Expression("u");
pic = new ParticleInitialConditionCount(initialCount, locationX, locationY, locationZ);
}
particleInitialConditions.add(pic);
//
// diffusion
//
Expression diffusion = new Expression(getMathSymbol(scs.getDiffusionParameter(), sm.getGeometryClass()));
Expression driftXExp = null;
if (scs.getVelocityXParameter().getExpression() != null) {
driftXExp = new Expression(getMathSymbol(scs.getVelocityXParameter(), sm.getGeometryClass()));
} else {
SpatialQuantity[] velX_quantities = scs.getVelocityQuantities(QuantityComponent.X);
if (velX_quantities.length > 0) {
int numRegions = simContext.getGeometry().getGeometrySurfaceDescription().getGeometricRegions(sm.getGeometryClass()).length;
if (velX_quantities.length == 1 && numRegions == 1) {
driftXExp = new Expression(getMathSymbol(velX_quantities[0], sm.getGeometryClass()));
} else {
throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
}
}
}
Expression driftYExp = null;
if (scs.getVelocityYParameter().getExpression() != null) {
driftYExp = new Expression(getMathSymbol(scs.getVelocityYParameter(), sm.getGeometryClass()));
} else {
SpatialQuantity[] velY_quantities = scs.getVelocityQuantities(QuantityComponent.Y);
if (velY_quantities.length > 0) {
int numRegions = simContext.getGeometry().getGeometrySurfaceDescription().getGeometricRegions(sm.getGeometryClass()).length;
if (velY_quantities.length == 1 && numRegions == 1) {
driftYExp = new Expression(getMathSymbol(velY_quantities[0], sm.getGeometryClass()));
} else {
throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
}
}
}
Expression driftZExp = null;
if (scs.getVelocityZParameter().getExpression() != null) {
driftZExp = new Expression(getMathSymbol(scs.getVelocityZParameter(), sm.getGeometryClass()));
} else {
SpatialQuantity[] velZ_quantities = scs.getVelocityQuantities(QuantityComponent.Z);
if (velZ_quantities.length > 0) {
int numRegions = simContext.getGeometry().getGeometrySurfaceDescription().getGeometricRegions(sm.getGeometryClass()).length;
if (velZ_quantities.length == 1 && numRegions == 1) {
driftZExp = new Expression(getMathSymbol(velZ_quantities[0], sm.getGeometryClass()));
} else {
throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
}
}
}
ParticleProperties particleProperties = new ParticleProperties(particleVariable, diffusion, driftXExp, driftYExp, driftZExp, particleInitialConditions);
GeometryClass myGC = sm.getGeometryClass();
if (myGC == null) {
throw new MappingException("Application '" + getSimulationContext().getName() + "'\nGeometry->StructureMapping->(" + sm.getStructure().getTypeName() + ")'" + sm.getStructure().getName() + "' must be mapped to geometry domain.\n(see 'Problems' tab)");
}
SubDomain subDomain = mathDesc.getSubDomain(myGC.getName());
subDomain.addParticleProperties(particleProperties);
}
}
for (ReactionStep reactionStep : reactionSteps) {
Kinetics kinetics = reactionStep.getKinetics();
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(reactionStep.getStructure());
GeometryClass reactionStepGeometryClass = sm.getGeometryClass();
SubDomain subdomain = mathDesc.getSubDomain(reactionStepGeometryClass.getName());
KineticsParameter reactionRateParameter = null;
if (kinetics instanceof LumpedKinetics) {
reactionRateParameter = ((LumpedKinetics) kinetics).getLumpedReactionRateParameter();
} else {
reactionRateParameter = ((DistributedKinetics) kinetics).getReactionRateParameter();
}
// macroscopic_irreversible/Microscopic_irreversible for bimolecular membrane reactions. They will NOT go through MassAction solver.
if (kinetics.getKineticsDescription().equals(KineticsDescription.Macroscopic_irreversible) || kinetics.getKineticsDescription().equals(KineticsDescription.Microscopic_irreversible)) {
Expression radiusExp = getIdentifierSubstitutions(reactionStep.getKinetics().getKineticsParameterFromRole(Kinetics.ROLE_Binding_Radius).getExpression(), modelUnitSystem.getBindingRadiusUnit(), reactionStepGeometryClass);
if (radiusExp != null) {
Expression expCopy = new Expression(radiusExp);
try {
MassActionSolver.substituteParameters(expCopy, true).evaluateConstant();
} catch (ExpressionException e) {
throw new MathException(VCellErrorMessages.getMassActionSolverMessage(reactionStep.getName(), "Problem in binding radius of " + reactionStep.getName() + ": '" + radiusExp.infix() + "', " + e.getMessage()));
}
} else {
throw new MathException(VCellErrorMessages.getMassActionSolverMessage(reactionStep.getName(), "Binding radius of " + reactionStep.getName() + " is null."));
}
List<ParticleVariable> reactantParticles = new ArrayList<ParticleVariable>();
List<ParticleVariable> productParticles = new ArrayList<ParticleVariable>();
List<Action> forwardActions = new ArrayList<Action>();
for (ReactionParticipant rp : reactionStep.getReactionParticipants()) {
SpeciesContext sc = rp.getSpeciesContext();
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
GeometryClass scGeometryClass = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass();
String varName = getMathSymbol(sc, scGeometryClass);
Variable var = mathDesc.getVariable(varName);
if (var instanceof ParticleVariable) {
ParticleVariable particle = (ParticleVariable) var;
if (rp instanceof Reactant) {
reactantParticles.add(particle);
if (!scs.isConstant() && !scs.isForceContinuous()) {
for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
if (radiusExp != null) {
forwardActions.add(Action.createDestroyAction(particle));
}
}
}
} else if (rp instanceof Product) {
productParticles.add(particle);
if (!scs.isConstant() && !scs.isForceContinuous()) {
for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
if (radiusExp != null) {
forwardActions.add(Action.createCreateAction(particle));
}
}
}
}
} else {
throw new MappingException("particle variable '" + varName + "' not found");
}
}
JumpProcessRateDefinition bindingRadius = new InteractionRadius(radiusExp);
// get jump process name
String jpName = TokenMangler.mangleToSName(reactionStep.getName());
// only for NFSim/Rules for now.
ProcessSymmetryFactor processSymmetryFactor = null;
if (forwardActions.size() > 0) {
ParticleJumpProcess forwardProcess = new ParticleJumpProcess(jpName, reactantParticles, bindingRadius, forwardActions, processSymmetryFactor);
subdomain.addParticleJumpProcess(forwardProcess);
}
} else // other type of reactions
{
/* 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;
// Using the MassActionFunction to write out the math description
MassActionSolver.MassActionFunction maFunc = null;
if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction) || kinetics.getKineticsDescription().equals(KineticsDescription.General) || kinetics.getKineticsDescription().equals(KineticsDescription.GeneralPermeability)) {
Expression rateExp = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate).getExpression();
Parameter forwardRateParameter = null;
Parameter reverseRateParameter = null;
if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction)) {
forwardRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KForward);
reverseRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KReverse);
} else if (kinetics.getKineticsDescription().equals(KineticsDescription.GeneralPermeability)) {
forwardRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_Permeability);
reverseRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_Permeability);
}
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();
}
}
}
if (maFunc != null) {
// 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
List<ParticleVariable> reactantParticles = new ArrayList<ParticleVariable>();
List<ParticleVariable> productParticles = new ArrayList<ParticleVariable>();
List<Action> forwardActions = new ArrayList<Action>();
List<Action> reverseActions = new ArrayList<Action>();
List<ReactionParticipant> reactants = maFunc.getReactants();
List<ReactionParticipant> products = maFunc.getProducts();
for (ReactionParticipant rp : reactants) {
SpeciesContext sc = rp.getSpeciesContext();
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
GeometryClass scGeometryClass = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass();
String varName = getMathSymbol(sc, scGeometryClass);
Variable var = mathDesc.getVariable(varName);
if (var instanceof ParticleVariable) {
ParticleVariable particle = (ParticleVariable) var;
reactantParticles.add(particle);
if (!scs.isConstant() && !scs.isForceContinuous()) {
for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
if (forwardRate != null) {
forwardActions.add(Action.createDestroyAction(particle));
}
if (reverseRate != null) {
reverseActions.add(Action.createCreateAction(particle));
}
}
}
} else {
throw new MappingException("particle variable '" + varName + "' not found");
}
}
for (ReactionParticipant rp : products) {
SpeciesContext sc = rp.getSpeciesContext();
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
GeometryClass scGeometryClass = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass();
String varName = getMathSymbol(sc, scGeometryClass);
Variable var = mathDesc.getVariable(varName);
if (var instanceof ParticleVariable) {
ParticleVariable particle = (ParticleVariable) var;
productParticles.add(particle);
if (!scs.isConstant() && !scs.isForceContinuous()) {
for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
if (forwardRate != null) {
forwardActions.add(Action.createCreateAction(particle));
}
if (reverseRate != null) {
reverseActions.add(Action.createDestroyAction(particle));
}
}
}
} else {
throw new MappingException("particle variable '" + varName + "' not found");
}
}
//
// There are two unit conversions required:
//
// 1) convert entire reaction rate from vcell reaction units to Smoldyn units (molecules/lengthunit^dim/timeunit)
// (where dim is 2 for membrane reactions and 3 for volume reactions)
//
// for forward rates:
// 2) convert each reactant from Smoldyn units (molecules/lengthunit^dim) to VCell units
// (where dim is 2 for membrane reactants and 3 for volume reactants)
//
// or
//
// for reverse rates:
// 2) convert each product from Smoldyn units (molecules/lengthunit^dim) to VCell units
// (where dim is 2 for membrane products and 3 for volume products)
//
RationalNumber reactionLocationDim = new RationalNumber(reactionStep.getStructure().getDimension());
VCUnitDefinition timeUnit = modelUnitSystem.getTimeUnit();
VCUnitDefinition smoldynReactionSizeUnit = modelUnitSystem.getLengthUnit().raiseTo(reactionLocationDim);
VCUnitDefinition smoldynSubstanceUnit = modelUnitSystem.getStochasticSubstanceUnit();
VCUnitDefinition smoldynReactionRateUnit = smoldynSubstanceUnit.divideBy(smoldynReactionSizeUnit).divideBy(timeUnit);
VCUnitDefinition vcellReactionRateUnit = reactionRateParameter.getUnitDefinition();
VCUnitDefinition reactionUnitFactor = smoldynReactionRateUnit.divideBy(vcellReactionRateUnit);
if (forwardRate != null) {
VCUnitDefinition smoldynReactantsUnit = modelUnitSystem.getInstance_DIMENSIONLESS();
// start with factor to translate entire reaction rate.
VCUnitDefinition forwardUnitFactor = reactionUnitFactor;
//
for (ReactionParticipant reactant : maFunc.getReactants()) {
VCUnitDefinition vcellReactantUnit = reactant.getSpeciesContext().getUnitDefinition();
boolean bForceContinuous = simContext.getReactionContext().getSpeciesContextSpec(reactant.getSpeciesContext()).isForceContinuous();
VCUnitDefinition smoldynReactantUnit = null;
if (bForceContinuous) {
// reactant is continuous (vcell units)
smoldynReactantUnit = reactant.getSpeciesContext().getUnitDefinition();
} else {
// reactant is a particle (smoldyn units)
RationalNumber reactantLocationDim = new RationalNumber(reactant.getStructure().getDimension());
VCUnitDefinition smoldynReactantSize = modelUnitSystem.getLengthUnit().raiseTo(reactantLocationDim);
smoldynReactantUnit = smoldynSubstanceUnit.divideBy(smoldynReactantSize);
}
// keep track of units of all reactants
smoldynReactantsUnit = smoldynReactantsUnit.multiplyBy(smoldynReactantUnit);
RationalNumber reactantStoichiometry = new RationalNumber(reactant.getStoichiometry());
VCUnitDefinition reactantUnitFactor = (vcellReactantUnit.divideBy(smoldynReactantUnit)).raiseTo(reactantStoichiometry);
// accumulate unit factors for all reactants
forwardUnitFactor = forwardUnitFactor.multiplyBy(reactantUnitFactor);
}
forwardRate = Expression.mult(forwardRate, getUnitFactor(forwardUnitFactor));
VCUnitDefinition smoldynExpectedForwardRateUnit = smoldynReactionRateUnit.divideBy(smoldynReactantsUnit);
// get probability
Expression exp = getIdentifierSubstitutions(forwardRate, smoldynExpectedForwardRateUnit, reactionStepGeometryClass).flatten();
JumpProcessRateDefinition partRateDef = new MacroscopicRateConstant(exp);
// create particle jump process
String jpName = TokenMangler.mangleToSName(reactionStep.getName());
// only for NFSim/Rules for now.
ProcessSymmetryFactor processSymmetryFactor = null;
if (forwardActions.size() > 0) {
ParticleJumpProcess forwardProcess = new ParticleJumpProcess(jpName, reactantParticles, partRateDef, forwardActions, processSymmetryFactor);
subdomain.addParticleJumpProcess(forwardProcess);
}
}
// end of forward rate not null
if (reverseRate != null) {
VCUnitDefinition smoldynProductsUnit = modelUnitSystem.getInstance_DIMENSIONLESS();
// start with factor to translate entire reaction rate.
VCUnitDefinition reverseUnitFactor = reactionUnitFactor;
//
for (ReactionParticipant product : maFunc.getProducts()) {
VCUnitDefinition vcellProductUnit = product.getSpeciesContext().getUnitDefinition();
boolean bForceContinuous = simContext.getReactionContext().getSpeciesContextSpec(product.getSpeciesContext()).isForceContinuous();
VCUnitDefinition smoldynProductUnit = null;
if (bForceContinuous) {
smoldynProductUnit = product.getSpeciesContext().getUnitDefinition();
} else {
RationalNumber productLocationDim = new RationalNumber(product.getStructure().getDimension());
VCUnitDefinition smoldynProductSize = modelUnitSystem.getLengthUnit().raiseTo(productLocationDim);
smoldynProductUnit = smoldynSubstanceUnit.divideBy(smoldynProductSize);
}
// keep track of units of all products
smoldynProductsUnit = smoldynProductsUnit.multiplyBy(smoldynProductUnit);
RationalNumber productStoichiometry = new RationalNumber(product.getStoichiometry());
VCUnitDefinition productUnitFactor = (vcellProductUnit.divideBy(smoldynProductUnit)).raiseTo(productStoichiometry);
// accumulate unit factors for all products
reverseUnitFactor = reverseUnitFactor.multiplyBy(productUnitFactor);
}
reverseRate = Expression.mult(reverseRate, getUnitFactor(reverseUnitFactor));
VCUnitDefinition smoldynExpectedReverseRateUnit = smoldynReactionRateUnit.divideBy(smoldynProductsUnit);
// get probability
Expression exp = getIdentifierSubstitutions(reverseRate, smoldynExpectedReverseRateUnit, reactionStepGeometryClass).flatten();
JumpProcessRateDefinition partProbRate = new MacroscopicRateConstant(exp);
// get jump process name
String jpName = TokenMangler.mangleToSName(reactionStep.getName() + "_reverse");
// only for NFSim/Rules for now.
ProcessSymmetryFactor processSymmetryFactor = null;
if (reverseActions.size() > 0) {
ParticleJumpProcess reverseProcess = new ParticleJumpProcess(jpName, productParticles, partProbRate, reverseActions, processSymmetryFactor);
subdomain.addParticleJumpProcess(reverseProcess);
}
}
// end of reverse rate not null
}
// end of maFunc not null
}
// end of reaction step for loop
}
//
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()) {
lg.warn(mathDesc.getVCML_database());
throw new MappingException("generated an invalid mathDescription: " + mathDesc.getWarning());
}
if (lg.isDebugEnabled()) {
System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string begin ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
System.out.println(mathDesc.getVCML());
System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string end ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
}
}
Also used :
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
LumpedKinetics(cbit.vcell.model.LumpedKinetics)
MathDescription(cbit.vcell.math.MathDescription)
ArrayList(java.util.ArrayList)
Product(cbit.vcell.model.Product)
SpeciesContext(cbit.vcell.model.SpeciesContext)
StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter)
Reactant(cbit.vcell.model.Reactant)
ExpressionException(cbit.vcell.parser.ExpressionException)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter)
SubVolume(cbit.vcell.geometry.SubVolume)
Vector(java.util.Vector)
SpatialQuantity(cbit.vcell.mapping.spatial.SpatialObject.SpatialQuantity)
SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter)
SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter)
JumpProcessRateDefinition(cbit.vcell.math.JumpProcessRateDefinition)
InteractionRadius(cbit.vcell.math.InteractionRadius)
ParticleJumpProcess(cbit.vcell.math.ParticleJumpProcess)
ModelParameter(cbit.vcell.model.Model.ModelParameter)
VCUnitDefinition(cbit.vcell.units.VCUnitDefinition)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
ParticleInitialCondition(cbit.vcell.math.ParticleProperties.ParticleInitialCondition)
ReactionStep(cbit.vcell.model.ReactionStep)
ParticleProperties(cbit.vcell.math.ParticleProperties)
Kinetics(cbit.vcell.model.Kinetics)
DistributedKinetics(cbit.vcell.model.DistributedKinetics)
LumpedKinetics(cbit.vcell.model.LumpedKinetics)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
Domain(cbit.vcell.math.Variable.Domain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
ReactionParticipant(cbit.vcell.model.ReactionParticipant)
GeometryClass(cbit.vcell.geometry.GeometryClass)
Action(cbit.vcell.math.Action)
VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable)
MembraneParticleVariable(cbit.vcell.math.MembraneParticleVariable)
ParticleVariable(cbit.vcell.math.ParticleVariable)
Variable(cbit.vcell.math.Variable)
SpeciesContextSpecProxyParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecProxyParameter)
SurfaceClass(cbit.vcell.geometry.SurfaceClass)
VariableHash(cbit.vcell.math.VariableHash)
VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable)
MembraneParticleVariable(cbit.vcell.math.MembraneParticleVariable)
ParticleVariable(cbit.vcell.math.ParticleVariable)
MacroscopicRateConstant(cbit.vcell.math.MacroscopicRateConstant)
Constant(cbit.vcell.math.Constant)
SymbolTableEntry(cbit.vcell.parser.SymbolTableEntry)
MacroscopicRateConstant(cbit.vcell.math.MacroscopicRateConstant)
RationalNumber(ucar.units_vcell.RationalNumber)
ModelUnitSystem(cbit.vcell.model.ModelUnitSystem)
Pair(org.vcell.util.Pair)
ParticleInitialConditionConcentration(cbit.vcell.math.ParticleProperties.ParticleInitialConditionConcentration)
ProcessSymmetryFactor(cbit.vcell.math.ParticleJumpProcess.ProcessSymmetryFactor)
Expression(cbit.vcell.parser.Expression)
VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable)
MathException(cbit.vcell.math.MathException)
Model(cbit.vcell.model.Model)
StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter)
Parameter(cbit.vcell.model.Parameter)
KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter)
SpeciesContextSpecProxyParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecProxyParameter)
ModelParameter(cbit.vcell.model.Model.ModelParameter)
SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter)
MassActionSolver(cbit.vcell.model.MassActionSolver)
ParticleInitialConditionCount(cbit.vcell.math.ParticleProperties.ParticleInitialConditionCount)
Example 8 with Reactant
use of cbit.vcell.model.Reactant in project vcell by virtualcell.
the class RulebasedTransformer method transform.
private void transform(SimulationContext originalSimContext, SimulationContext transformedSimulationContext, ArrayList<ModelEntityMapping> entityMappings, MathMappingCallback mathMappingCallback) throws PropertyVetoException {
Model newModel = transformedSimulationContext.getModel();
Model originalModel = originalSimContext.getModel();
ModelEntityMapping em = null;
// list of rules created from the reactions; we apply the symmetry factor computed by bionetgen only to these
Set<ReactionRule> fromReactions = new HashSet<>();
for (SpeciesContext newSpeciesContext : newModel.getSpeciesContexts()) {
final SpeciesContext originalSpeciesContext = originalModel.getSpeciesContext(newSpeciesContext.getName());
// map new and old species contexts
em = new ModelEntityMapping(originalSpeciesContext, newSpeciesContext);
entityMappings.add(em);
if (newSpeciesContext.hasSpeciesPattern()) {
// it's perfect already and can't be improved
continue;
}
try {
MolecularType newmt = newModel.getRbmModelContainer().createMolecularType();
newModel.getRbmModelContainer().addMolecularType(newmt, false);
MolecularTypePattern newmtp_sc = new MolecularTypePattern(newmt);
SpeciesPattern newsp_sc = new SpeciesPattern();
newsp_sc.addMolecularTypePattern(newmtp_sc);
newSpeciesContext.setSpeciesPattern(newsp_sc);
RbmObservable newo = new RbmObservable(newModel, "O0_" + newmt.getName() + "_tot", newSpeciesContext.getStructure(), RbmObservable.ObservableType.Molecules);
MolecularTypePattern newmtp_ob = new MolecularTypePattern(newmt);
SpeciesPattern newsp_ob = new SpeciesPattern();
newsp_ob.addMolecularTypePattern(newmtp_ob);
newo.addSpeciesPattern(newsp_ob);
newModel.getRbmModelContainer().addObservable(newo);
// map new observable to old species context
em = new ModelEntityMapping(originalSpeciesContext, newo);
entityMappings.add(em);
} catch (ModelException e) {
e.printStackTrace();
throw new RuntimeException("unable to transform species context: " + e.getMessage());
} catch (PropertyVetoException e) {
// TODO Auto-generated catch block
e.printStackTrace();
}
}
ReactionSpec[] reactionSpecs = transformedSimulationContext.getReactionContext().getReactionSpecs();
for (ReactionSpec reactionSpec : reactionSpecs) {
if (reactionSpec.isExcluded()) {
// we create rules only from those reactions which are not excluded
continue;
}
ReactionStep rs = reactionSpec.getReactionStep();
String name = rs.getName();
String mangled = TokenMangler.fixTokenStrict(name);
mangled = newModel.getReactionName(mangled);
Kinetics k = rs.getKinetics();
if (!(k instanceof MassActionKinetics)) {
throw new RuntimeException("Only Mass Action Kinetics supported at this time, reaction \"" + rs.getName() + "\" uses kinetic law type \"" + rs.getKinetics().getName() + "\"");
}
boolean bReversible = rs.isReversible();
ReactionRule rr = new ReactionRule(newModel, mangled, rs.getStructure(), bReversible);
fromReactions.add(rr);
MassActionKinetics massActionKinetics = (MassActionKinetics) k;
List<Reactant> rList = rs.getReactants();
List<Product> pList = rs.getProducts();
// counting the stoichiometry - 2A+B means 3 reactants
int numReactants = 0;
for (Reactant r : rList) {
numReactants += r.getStoichiometry();
if (numReactants > 2) {
String message = "NFSim doesn't support more than 2 reactants within a reaction: " + name;
throw new RuntimeException(message);
}
}
int numProducts = 0;
for (Product p : pList) {
numProducts += p.getStoichiometry();
if (bReversible && numProducts > 2) {
String message = "NFSim doesn't support more than 2 products within a reversible reaction: " + name;
throw new RuntimeException(message);
}
}
RateLawType rateLawType = RateLawType.MassAction;
RbmKineticLaw kineticLaw = new RbmKineticLaw(rr, rateLawType);
try {
String forwardRateName = massActionKinetics.getForwardRateParameter().getName();
Expression forwardRateExp = massActionKinetics.getForwardRateParameter().getExpression();
String reverseRateName = massActionKinetics.getReverseRateParameter().getName();
Expression reverseRateExp = massActionKinetics.getReverseRateParameter().getExpression();
LocalParameter fR = kineticLaw.getLocalParameter(RbmKineticLawParameterType.MassActionForwardRate);
fR.setName(forwardRateName);
LocalParameter rR = kineticLaw.getLocalParameter(RbmKineticLawParameterType.MassActionReverseRate);
rR.setName(reverseRateName);
if (rs.hasReactant()) {
kineticLaw.setParameterValue(fR, forwardRateExp, true);
}
if (rs.hasProduct()) {
kineticLaw.setParameterValue(rR, reverseRateExp, true);
}
//
for (KineticsParameter reaction_p : massActionKinetics.getKineticsParameters()) {
if (reaction_p.getRole() == Kinetics.ROLE_UserDefined) {
LocalParameter rule_p = kineticLaw.getLocalParameter(reaction_p.getName());
if (rule_p == null) {
//
// after lazy parameter creation we didn't find a user-defined rule parameter with this same name.
//
// there must be a global symbol with the same name, that the local reaction parameter has overridden.
//
ParameterContext.LocalProxyParameter rule_proxy_parameter = null;
for (ProxyParameter proxyParameter : kineticLaw.getProxyParameters()) {
if (proxyParameter.getName().equals(reaction_p.getName())) {
rule_proxy_parameter = (LocalProxyParameter) proxyParameter;
}
}
if (rule_proxy_parameter != null) {
// we want to convert to local
boolean bConvertToGlobal = false;
kineticLaw.convertParameterType(rule_proxy_parameter, bConvertToGlobal);
} else {
// could find neither local parameter nor proxy parameter
throw new RuntimeException("user defined parameter " + reaction_p.getName() + " from reaction " + rs.getName() + " didn't map to a reactionRule parameter");
}
} else if (rule_p.getRole() == RbmKineticLawParameterType.UserDefined) {
kineticLaw.setParameterValue(rule_p, reaction_p.getExpression(), true);
rule_p.setUnitDefinition(reaction_p.getUnitDefinition());
} else {
throw new RuntimeException("user defined parameter " + reaction_p.getName() + " from reaction " + rs.getName() + " mapped to a reactionRule parameter with unexpected role " + rule_p.getRole().getDescription());
}
}
}
} catch (ExpressionException e) {
e.printStackTrace();
throw new RuntimeException("Problem attempting to set RbmKineticLaw expression: " + e.getMessage());
}
rr.setKineticLaw(kineticLaw);
KineticsParameter[] kpList = k.getKineticsParameters();
ModelParameter[] mpList = rs.getModel().getModelParameters();
ModelParameter mp = rs.getModel().getModelParameter(kpList[0].getName());
ReactionParticipant[] rpList = rs.getReactionParticipants();
for (ReactionParticipant p : rpList) {
if (p instanceof Reactant) {
int stoichiometry = p.getStoichiometry();
for (int i = 0; i < stoichiometry; i++) {
SpeciesPattern speciesPattern = new SpeciesPattern(rs.getModel(), p.getSpeciesContext().getSpeciesPattern());
ReactantPattern reactantPattern = new ReactantPattern(speciesPattern, p.getStructure());
rr.addReactant(reactantPattern);
}
} else if (p instanceof Product) {
int stoichiometry = p.getStoichiometry();
for (int i = 0; i < stoichiometry; i++) {
SpeciesPattern speciesPattern = new SpeciesPattern(rs.getModel(), p.getSpeciesContext().getSpeciesPattern());
ProductPattern productPattern = new ProductPattern(speciesPattern, p.getStructure());
rr.addProduct(productPattern);
}
}
}
// commented code below is probably obsolete, we verify (above) in the reaction the number of participants,
// no need to do it again in the corresponding rule
// if(rr.getReactantPatterns().size() > 2) {
// String message = "NFSim doesn't support more than 2 reactants within a reaction: " + name;
// throw new RuntimeException(message);
// }
// if(rr.getProductPatterns().size() > 2) {
// String message = "NFSim doesn't support more than 2 products within a reaction: " + name;
// throw new RuntimeException(message);
// }
newModel.removeReactionStep(rs);
newModel.getRbmModelContainer().addReactionRule(rr);
}
for (ReactionRuleSpec rrs : transformedSimulationContext.getReactionContext().getReactionRuleSpecs()) {
if (rrs == null) {
continue;
}
ReactionRule rr = rrs.getReactionRule();
if (rrs.isExcluded()) {
// delete those rules which are disabled (excluded) in the Specifications / Reaction table
newModel.getRbmModelContainer().removeReactionRule(rr);
continue;
}
}
// now that we generated the rules we can delete the reaction steps they're coming from
for (ReactionStep rs : newModel.getReactionSteps()) {
newModel.removeReactionStep(rs);
}
try {
// we invoke bngl just for the purpose of generating the xml file, which we'll then use to extract the symmetry factor
generateNetwork(transformedSimulationContext, fromReactions, mathMappingCallback);
} catch (ClassNotFoundException | IOException e) {
// TODO Auto-generated catch block
e.printStackTrace();
}
System.out.println("Finished RuleBased Transformer.");
}
Also used :
Product(cbit.vcell.model.Product)
SpeciesContext(cbit.vcell.model.SpeciesContext)
Reactant(cbit.vcell.model.Reactant)
LocalProxyParameter(cbit.vcell.mapping.ParameterContext.LocalProxyParameter)
SpeciesPattern(org.vcell.model.rbm.SpeciesPattern)
ExpressionException(cbit.vcell.parser.ExpressionException)
KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter)
RateLawType(cbit.vcell.model.RbmKineticLaw.RateLawType)
HashSet(java.util.HashSet)
ReactantPattern(cbit.vcell.model.ReactantPattern)
ReactionRule(cbit.vcell.model.ReactionRule)
ModelException(cbit.vcell.model.ModelException)
ProductPattern(cbit.vcell.model.ProductPattern)
RbmObservable(cbit.vcell.model.RbmObservable)
RbmKineticLaw(cbit.vcell.model.RbmKineticLaw)
IOException(java.io.IOException)
MolecularType(org.vcell.model.rbm.MolecularType)
PropertyVetoException(java.beans.PropertyVetoException)
LocalParameter(cbit.vcell.mapping.ParameterContext.LocalParameter)
ModelParameter(cbit.vcell.model.Model.ModelParameter)
ProxyParameter(cbit.vcell.model.ProxyParameter)
LocalProxyParameter(cbit.vcell.mapping.ParameterContext.LocalProxyParameter)
Expression(cbit.vcell.parser.Expression)
ReactionStep(cbit.vcell.model.ReactionStep)
Model(cbit.vcell.model.Model)
MassActionKinetics(cbit.vcell.model.MassActionKinetics)
Kinetics(cbit.vcell.model.Kinetics)
MassActionKinetics(cbit.vcell.model.MassActionKinetics)
MolecularTypePattern(org.vcell.model.rbm.MolecularTypePattern)
ReactionParticipant(cbit.vcell.model.ReactionParticipant)
Example 9 with Reactant
use of cbit.vcell.model.Reactant in project vcell by virtualcell.
the class NetworkTransformer method transform.
private void transform(SimulationContext simContext, SimulationContext transformedSimulationContext, ArrayList<ModelEntityMapping> entityMappings, MathMappingCallback mathMappingCallback, NetworkGenerationRequirements networkGenerationRequirements) {
String msg = "Generating network: flattening...";
mathMappingCallback.setMessage(msg);
TaskCallbackMessage tcm = new TaskCallbackMessage(TaskCallbackStatus.Clean, "");
simContext.appendToConsole(tcm);
tcm = new TaskCallbackMessage(TaskCallbackStatus.TaskStart, msg);
simContext.appendToConsole(tcm);
long startTime = System.currentTimeMillis();
System.out.println("Convert to bngl, execute BNG, retrieve the results.");
try {
BNGOutputSpec outputSpec = generateNetwork(simContext, mathMappingCallback, networkGenerationRequirements);
if (mathMappingCallback.isInterrupted()) {
msg = "Canceled by user.";
tcm = new TaskCallbackMessage(TaskCallbackStatus.Error, msg);
simContext.appendToConsole(tcm);
throw new UserCancelException(msg);
}
long endTime = System.currentTimeMillis();
long elapsedTime = endTime - startTime;
System.out.println(" " + elapsedTime + " milliseconds");
Model model = transformedSimulationContext.getModel();
ReactionContext reactionContext = transformedSimulationContext.getReactionContext();
// ---- Parameters -----------------------------------------------------------------------------------------------
startTime = System.currentTimeMillis();
for (int i = 0; i < outputSpec.getBNGParams().length; i++) {
BNGParameter p = outputSpec.getBNGParams()[i];
// System.out.println(i+1 + ":\t\t"+ p.toString());
if (model.getRbmModelContainer().getParameter(p.getName()) != null) {
// if it's already there we don't try to add it again; this should be true for all of them!
continue;
}
String s = p.getName();
if (NetworkConstraints.SPECIES_LIMIT_PARAMETER.equals(s)) {
System.out.println("found NetworkConstraints seciesLimit parameter.");
continue;
}
if (NetworkConstraints.REACTIONS_LIMIT_PARAMETER.equals(s)) {
System.out.println("found NetworkConstraints reactionsLimit parameter.");
continue;
}
FakeSeedSpeciesInitialConditionsParameter fakeICParam = FakeSeedSpeciesInitialConditionsParameter.fromString(s);
if (speciesEquivalenceMap.containsKey(fakeICParam)) {
// we get rid of the fake parameters we use as keys
continue;
}
FakeReactionRuleRateParameter fakeKineticParam = FakeReactionRuleRateParameter.fromString(s);
if (fakeKineticParam != null) {
System.out.println("found fakeKineticParam " + fakeKineticParam.fakeParameterName);
// we get rid of the fake parameters we use as keys
continue;
}
throw new RuntimeException("unexpected parameter " + p.getName() + " in internal BNG processing");
// Expression exp = new Expression(p.getValue());
// exp.bindExpression(model.getRbmModelContainer().getSymbolTable());
// model.getRbmModelContainer().addParameter(p.getName(), exp, model.getUnitSystem().getInstance_TBD());
}
endTime = System.currentTimeMillis();
elapsedTime = endTime - startTime;
msg = "Adding " + outputSpec.getBNGParams().length + " parameters to model, " + elapsedTime + " ms";
System.out.println(msg);
// ---- Species ------------------------------------------------------------------------------------------------------------
mathMappingCallback.setMessage("generating network: adding species...");
mathMappingCallback.setProgressFraction(progressFractionQuota / 4.0f);
startTime = System.currentTimeMillis();
System.out.println("\nSpecies :");
// the reactions will need this map to recover the names of species knowing only the networkFileIndex
HashMap<Integer, String> speciesMap = new HashMap<Integer, String>();
LinkedHashMap<String, Species> sMap = new LinkedHashMap<String, Species>();
LinkedHashMap<String, SpeciesContext> scMap = new LinkedHashMap<String, SpeciesContext>();
LinkedHashMap<String, BNGSpecies> crossMap = new LinkedHashMap<String, BNGSpecies>();
List<SpeciesContext> noMapForThese = new ArrayList<SpeciesContext>();
// final int decimalTickCount = Math.max(outputSpec.getBNGSpecies().length/10, 1);
for (int i = 0; i < outputSpec.getBNGSpecies().length; i++) {
BNGSpecies s = outputSpec.getBNGSpecies()[i];
// System.out.println(i+1 + ":\t\t"+ s.toString());
String key = s.getConcentration().infix();
FakeSeedSpeciesInitialConditionsParameter fakeParam = FakeSeedSpeciesInitialConditionsParameter.fromString(key);
if (fakeParam != null) {
Pair<SpeciesContext, Expression> value = speciesEquivalenceMap.get(fakeParam);
// the species context of the original model
SpeciesContext originalsc = value.one;
Expression initial = value.two;
// replace the fake initial condition with the real one
s.setConcentration(initial);
// we'll have to find the species context from the cloned model which correspond to the original species
SpeciesContext sc = model.getSpeciesContext(originalsc.getName());
// System.out.println(sc.getName() + ", " + sc.getSpecies().getCommonName() + " ...is one of the original seed species.");
// existing name
speciesMap.put(s.getNetworkFileIndex(), sc.getName());
sMap.put(sc.getName(), sc.getSpecies());
scMap.put(sc.getName(), sc);
crossMap.put(sc.getName(), s);
noMapForThese.add(sc);
continue;
}
// all these species are new!
// generate unique name for the species
int count = 0;
String speciesName = null;
String nameRoot = "s";
String speciesPatternNameString = s.extractName();
while (true) {
speciesName = nameRoot + count;
if (Model.isNameUnused(speciesName, model) && !sMap.containsKey(speciesName) && !scMap.containsKey(speciesName)) {
break;
}
count++;
}
// newly created name
speciesMap.put(s.getNetworkFileIndex(), speciesName);
SpeciesContext speciesContext;
if (s.hasCompartment()) {
String speciesPatternCompartmentString = s.extractCompartment();
speciesContext = new SpeciesContext(new Species(speciesName, s.getName()), model.getStructure(speciesPatternCompartmentString), null);
} else {
speciesContext = new SpeciesContext(new Species(speciesName, s.getName()), model.getStructure(0), null);
}
speciesContext.setName(speciesName);
try {
if (speciesPatternNameString != null) {
SpeciesPattern sp = RbmUtils.parseSpeciesPattern(speciesPatternNameString, model);
speciesContext.setSpeciesPattern(sp);
}
} catch (ParseException e) {
e.printStackTrace();
throw new RuntimeException("Bad format for species pattern string: " + e.getMessage());
}
// speciesContext.setSpeciesPatternString(speciesPatternString);
// model.addSpecies(speciesContext.getSpecies());
// model.addSpeciesContext(speciesContext);
sMap.put(speciesName, speciesContext.getSpecies());
scMap.put(speciesName, speciesContext);
crossMap.put(speciesName, s);
// }
if (mathMappingCallback.isInterrupted()) {
msg = "Canceled by user.";
tcm = new TaskCallbackMessage(TaskCallbackStatus.Error, msg);
simContext.appendToConsole(tcm);
throw new UserCancelException(msg);
}
// if(i%50 == 0) {
// System.out.println(i+"");
// }
// if(i%decimalTickCount == 0) {
// int multiplier = i/decimalTickCount;
// float progress = progressFractionQuota/4.0f + progressFractionQuotaSpecies*multiplier;
// mathMappingCallback.setProgressFraction(progress);
// }
}
for (SpeciesContext sc1 : model.getSpeciesContexts()) {
boolean found = false;
for (Map.Entry<String, SpeciesContext> entry : scMap.entrySet()) {
SpeciesContext sc2 = entry.getValue();
if (sc1.getName().equals(sc2.getName())) {
found = true;
// System.out.println("found species context " + sc1.getName() + " of species " + sc1.getSpecies().getCommonName() + " // " + sc2.getSpecies().getCommonName());
break;
}
}
if (found == false) {
// we add to the map the species context and the species which exist in the model but which are not in the map yet
// the only ones in this situation should be plain species which were not given to bngl for flattening (they are flat already)
// System.out.println("species context " + sc1.getName() + " not found in the map. Adding it.");
scMap.put(sc1.getName(), sc1);
sMap.put(sc1.getName(), sc1.getSpecies());
noMapForThese.add(sc1);
}
}
for (Species s1 : model.getSpecies()) {
boolean found = false;
for (Map.Entry<String, Species> entry : sMap.entrySet()) {
Species s2 = entry.getValue();
if (s1.getCommonName().equals(s2.getCommonName())) {
found = true;
// System.out.println("found species " + s1.getCommonName());
break;
}
}
if (found == false) {
System.err.println("species " + s1.getCommonName() + " not found in the map!");
}
}
SpeciesContext[] sca = new SpeciesContext[scMap.size()];
scMap.values().toArray(sca);
Species[] sa = new HashSet<Species>(sMap.values()).toArray(new Species[0]);
model.setSpecies(sa);
model.setSpeciesContexts(sca);
boolean isSpatial = transformedSimulationContext.getGeometry().getDimension() > 0;
for (SpeciesContext sc : sca) {
if (noMapForThese.contains(sc)) {
continue;
}
SpeciesContextSpec scs = reactionContext.getSpeciesContextSpec(sc);
Parameter param = scs.getParameter(SpeciesContextSpec.ROLE_InitialConcentration);
BNGSpecies s = crossMap.get(sc.getName());
param.setExpression(s.getConcentration());
SpeciesContext origSpeciesContext = simContext.getModel().getSpeciesContext(s.getName());
if (origSpeciesContext != null) {
ModelEntityMapping em = new ModelEntityMapping(origSpeciesContext, sc);
entityMappings.add(em);
} else {
ModelEntityMapping em = new ModelEntityMapping(new GeneratedSpeciesSymbolTableEntry(sc), sc);
if (isSpatial) {
scs.initializeForSpatial();
}
entityMappings.add(em);
}
}
// for(SpeciesContext sc : sca) { // clean all the species patterns from the flattened species, we have no sp now
// sc.setSpeciesPattern(null);
// }
endTime = System.currentTimeMillis();
elapsedTime = endTime - startTime;
msg = "Adding " + outputSpec.getBNGSpecies().length + " species to model, " + elapsedTime + " ms";
System.out.println(msg);
// ---- Reactions -----------------------------------------------------------------------------------------------------
mathMappingCallback.setMessage("generating network: adding reactions...");
mathMappingCallback.setProgressFraction(progressFractionQuota / 4.0f * 3.0f);
startTime = System.currentTimeMillis();
System.out.println("\nReactions :");
Map<String, HashSet<String>> ruleKeyMap = new HashMap<String, HashSet<String>>();
Map<String, BNGReaction> directBNGReactionsMap = new HashMap<String, BNGReaction>();
Map<String, BNGReaction> reverseBNGReactionsMap = new HashMap<String, BNGReaction>();
for (int i = 0; i < outputSpec.getBNGReactions().length; i++) {
BNGReaction r = outputSpec.getBNGReactions()[i];
if (!r.isRuleReversed()) {
// direct
directBNGReactionsMap.put(r.getKey(), r);
} else {
reverseBNGReactionsMap.put(r.getKey(), r);
}
//
// for each rule name, store set of keySets (number of unique keysets are number of generated reactions from this ruleName).
//
HashSet<String> keySet = ruleKeyMap.get(r.getRuleName());
if (keySet == null) {
keySet = new HashSet<String>();
ruleKeyMap.put(r.getRuleName(), keySet);
}
keySet.add(r.getKey());
}
Map<String, ReactionStep> reactionStepMap = new HashMap<String, ReactionStep>();
for (int i = 0; i < outputSpec.getBNGReactions().length; i++) {
BNGReaction bngReaction = outputSpec.getBNGReactions()[i];
// System.out.println(i+1 + ":\t\t"+ r.writeReaction());
String baseName = bngReaction.getRuleName();
// which may never happen since we encountered this problem only one time in many years
if (baseName.contains(",") && (baseName.length() > 192)) {
int pos = baseName.indexOf(",");
baseName = baseName.substring(0, pos);
}
// System.out.println(i + ": " + baseName);
String reactionName = null;
HashSet<String> keySetsForThisRule = ruleKeyMap.get(bngReaction.getRuleName());
if (keySetsForThisRule.size() == 1 && model.getReactionStep(bngReaction.getRuleName()) == null && !reactionStepMap.containsKey(bngReaction.getRuleName())) {
// we can reuse the reaction rule labels
reactionName = baseName;
} else {
reactionName = baseName + "_0";
while (true) {
if (model.getReactionStep(reactionName) == null && !reactionStepMap.containsKey(reactionName)) {
// we can reuse the reaction rule labels
break;
}
reactionName = TokenMangler.getNextEnumeratedToken(reactionName);
}
}
//
if (directBNGReactionsMap.containsValue(bngReaction)) {
BNGReaction forwardBNGReaction = bngReaction;
BNGReaction reverseBNGReaction = reverseBNGReactionsMap.get(bngReaction.getKey());
String name = forwardBNGReaction.getRuleName();
if (name.endsWith(ReactionRule.DirectHalf)) {
name = name.substring(0, name.indexOf(ReactionRule.DirectHalf));
}
if (name.endsWith(ReactionRule.InverseHalf)) {
name = name.substring(0, name.indexOf(ReactionRule.InverseHalf));
}
ReactionRule rr = model.getRbmModelContainer().getReactionRule(name);
Structure structure = rr.getStructure();
boolean bReversible = reverseBNGReaction != null;
SimpleReaction sr = new SimpleReaction(model, structure, reactionName, bReversible);
for (int j = 0; j < forwardBNGReaction.getReactants().length; j++) {
BNGSpecies s = forwardBNGReaction.getReactants()[j];
String scName = speciesMap.get(s.getNetworkFileIndex());
SpeciesContext sc = model.getSpeciesContext(scName);
Reactant reactant = sr.getReactant(scName);
if (reactant == null) {
int stoichiometry = 1;
sr.addReactant(sc, stoichiometry);
} else {
int stoichiometry = reactant.getStoichiometry();
stoichiometry += 1;
reactant.setStoichiometry(stoichiometry);
}
}
for (int j = 0; j < forwardBNGReaction.getProducts().length; j++) {
BNGSpecies s = forwardBNGReaction.getProducts()[j];
String scName = speciesMap.get(s.getNetworkFileIndex());
SpeciesContext sc = model.getSpeciesContext(scName);
Product product = sr.getProduct(scName);
if (product == null) {
int stoichiometry = 1;
sr.addProduct(sc, stoichiometry);
} else {
int stoichiometry = product.getStoichiometry();
stoichiometry += 1;
product.setStoichiometry(stoichiometry);
}
}
if (!bngReaction.isMichaelisMenten()) {
// MassAction
MassActionKinetics targetKinetics = new MassActionKinetics(sr);
sr.setKinetics(targetKinetics);
KineticsParameter kforward = targetKinetics.getForwardRateParameter();
KineticsParameter kreverse = targetKinetics.getReverseRateParameter();
String kforwardNewName = rr.getKineticLaw().getLocalParameter(RbmKineticLawParameterType.MassActionForwardRate).getName();
if (!kforward.getName().equals(kforwardNewName)) {
targetKinetics.renameParameter(kforward.getName(), kforwardNewName);
kforward = targetKinetics.getForwardRateParameter();
}
final String kreverseNewName = rr.getKineticLaw().getLocalParameter(RbmKineticLawParameterType.MassActionReverseRate).getName();
if (!kreverse.getName().equals(kreverseNewName)) {
targetKinetics.renameParameter(kreverse.getName(), kreverseNewName);
kreverse = targetKinetics.getReverseRateParameter();
}
applyKineticsExpressions(forwardBNGReaction, kforward, targetKinetics);
if (reverseBNGReaction != null) {
applyKineticsExpressions(reverseBNGReaction, kreverse, targetKinetics);
}
} else {
// MichaelisMenten
HMM_IRRKinetics targetKinetics = new HMM_IRRKinetics(sr);
sr.setKinetics(targetKinetics);
KineticsParameter vmax = targetKinetics.getVmaxParameter();
KineticsParameter km = targetKinetics.getKmParameter();
String vmaxNewName = rr.getKineticLaw().getLocalParameter(RbmKineticLawParameterType.MichaelisMentenVmax).getName();
if (!vmax.getName().equals(vmaxNewName)) {
targetKinetics.renameParameter(vmax.getName(), vmaxNewName);
vmax = targetKinetics.getVmaxParameter();
}
applyKineticsExpressions(forwardBNGReaction, vmax, targetKinetics);
final String kmNewName = rr.getKineticLaw().getLocalParameter(RbmKineticLawParameterType.MichaelisMentenKm).getName();
if (!km.getName().equals(kmNewName)) {
targetKinetics.renameParameter(km.getName(), kmNewName);
km = targetKinetics.getKmParameter();
}
applyKineticsExpressions(forwardBNGReaction, km, targetKinetics);
}
reactionStepMap.put(reactionName, sr);
} else if (reverseBNGReactionsMap.containsValue(bngReaction) && !directBNGReactionsMap.containsKey(bngReaction.getKey())) {
// reverse only (must be irreversible, cannot be Michaelis-Menten)
BNGReaction reverseBNGReaction = reverseBNGReactionsMap.get(bngReaction.getKey());
ReactionRule rr = model.getRbmModelContainer().getReactionRule(reverseBNGReaction.extractRuleName());
Structure structure = rr.getStructure();
boolean bReversible = false;
SimpleReaction sr = new SimpleReaction(model, structure, reactionName, bReversible);
for (int j = 0; j < reverseBNGReaction.getReactants().length; j++) {
BNGSpecies s = reverseBNGReaction.getReactants()[j];
String scName = speciesMap.get(s.getNetworkFileIndex());
SpeciesContext sc = model.getSpeciesContext(scName);
Reactant reactant = sr.getReactant(scName);
if (reactant == null) {
int stoichiometry = 1;
sr.addReactant(sc, stoichiometry);
} else {
int stoichiometry = reactant.getStoichiometry();
stoichiometry += 1;
reactant.setStoichiometry(stoichiometry);
}
}
for (int j = 0; j < reverseBNGReaction.getProducts().length; j++) {
BNGSpecies s = reverseBNGReaction.getProducts()[j];
String scName = speciesMap.get(s.getNetworkFileIndex());
SpeciesContext sc = model.getSpeciesContext(scName);
Product product = sr.getProduct(scName);
if (product == null) {
int stoichiometry = 1;
sr.addProduct(sc, stoichiometry);
} else {
int stoichiometry = product.getStoichiometry();
stoichiometry += 1;
product.setStoichiometry(stoichiometry);
}
}
MassActionKinetics k = new MassActionKinetics(sr);
sr.setKinetics(k);
KineticsParameter kforward = k.getForwardRateParameter();
KineticsParameter kreverse = k.getReverseRateParameter();
String kforwardNewName = rr.getKineticLaw().getLocalParameter(RbmKineticLawParameterType.MassActionForwardRate).getName();
if (!kforward.getName().equals(kforwardNewName)) {
k.renameParameter(kforward.getName(), kforwardNewName);
kforward = k.getForwardRateParameter();
}
final String kreverseNewName = rr.getKineticLaw().getLocalParameter(RbmKineticLawParameterType.MassActionReverseRate).getName();
if (!kreverse.getName().equals(kreverseNewName)) {
k.renameParameter(kreverse.getName(), kreverseNewName);
kreverse = k.getReverseRateParameter();
}
applyKineticsExpressions(reverseBNGReaction, kforward, k);
// String fieldParameterName = kforward.getName();
// fieldParameterName += "_" + r.getRuleName();
// kforward.setName(fieldParameterName);
reactionStepMap.put(reactionName, sr);
}
}
for (ReactionStep rs : model.getReactionSteps()) {
reactionStepMap.put(rs.getName(), rs);
}
ReactionStep[] reactionSteps = new ReactionStep[reactionStepMap.size()];
reactionStepMap.values().toArray(reactionSteps);
model.setReactionSteps(reactionSteps);
if (mathMappingCallback.isInterrupted()) {
msg = "Canceled by user.";
tcm = new TaskCallbackMessage(TaskCallbackStatus.Error, msg);
simContext.appendToConsole(tcm);
throw new UserCancelException(msg);
}
endTime = System.currentTimeMillis();
elapsedTime = endTime - startTime;
msg = "Adding " + outputSpec.getBNGReactions().length + " reactions to model, " + elapsedTime + " ms";
System.out.println(msg);
// clean all the reaction rules
model.getRbmModelContainer().getReactionRuleList().clear();
// ---- Observables -------------------------------------------------------------------------------------------------
mathMappingCallback.setMessage("generating network: adding observables...");
mathMappingCallback.setProgressFraction(progressFractionQuota / 8.0f * 7.0f);
startTime = System.currentTimeMillis();
System.out.println("\nObservables :");
RbmModelContainer rbmmc = model.getRbmModelContainer();
for (int i = 0; i < outputSpec.getObservableGroups().length; i++) {
ObservableGroup o = outputSpec.getObservableGroups()[i];
if (rbmmc.getParameter(o.getObservableGroupName()) != null) {
System.out.println(" ...already exists.");
// if it's already there we don't try to add it again; this should be true for all of them!
continue;
}
ArrayList<Expression> terms = new ArrayList<Expression>();
for (int j = 0; j < o.getListofSpecies().length; j++) {
Expression term = Expression.mult(new Expression(o.getSpeciesMultiplicity()[j]), new Expression(speciesMap.get(o.getListofSpecies()[j].getNetworkFileIndex())));
terms.add(term);
}
Expression exp = Expression.add(terms.toArray(new Expression[terms.size()])).flatten();
exp.bindExpression(rbmmc.getSymbolTable());
RbmObservable originalObservable = rbmmc.getObservable(o.getObservableGroupName());
VCUnitDefinition observableUnitDefinition = originalObservable.getUnitDefinition();
rbmmc.removeObservable(originalObservable);
Parameter newParameter = rbmmc.addParameter(o.getObservableGroupName(), exp, observableUnitDefinition);
RbmObservable origObservable = simContext.getModel().getRbmModelContainer().getObservable(o.getObservableGroupName());
ModelEntityMapping em = new ModelEntityMapping(origObservable, newParameter);
entityMappings.add(em);
}
if (mathMappingCallback.isInterrupted()) {
msg = "Canceled by user.";
tcm = new TaskCallbackMessage(TaskCallbackStatus.Error, msg);
simContext.appendToConsole(tcm);
throw new UserCancelException(msg);
}
endTime = System.currentTimeMillis();
elapsedTime = endTime - startTime;
msg = "Adding " + outputSpec.getObservableGroups().length + " observables to model, " + elapsedTime + " ms";
System.out.println(msg);
} catch (PropertyVetoException ex) {
ex.printStackTrace(System.out);
throw new RuntimeException(ex.getMessage());
} catch (ExpressionBindingException ex) {
ex.printStackTrace(System.out);
throw new RuntimeException(ex.getMessage());
} catch (ModelException ex) {
ex.printStackTrace(System.out);
throw new RuntimeException(ex.getMessage());
} catch (ExpressionException ex) {
ex.printStackTrace(System.out);
throw new RuntimeException(ex.getMessage());
} catch (ClassNotFoundException ex) {
throw new RuntimeException(ex.getMessage());
} catch (IOException ex) {
throw new RuntimeException(ex.getMessage());
}
System.out.println("Done transforming");
msg = "Generating math...";
System.out.println(msg);
mathMappingCallback.setMessage(msg);
mathMappingCallback.setProgressFraction(progressFractionQuota);
}
Also used :
HashMap(java.util.HashMap)
LinkedHashMap(java.util.LinkedHashMap)
HMM_IRRKinetics(cbit.vcell.model.HMM_IRRKinetics)
UserCancelException(org.vcell.util.UserCancelException)
ArrayList(java.util.ArrayList)
Product(cbit.vcell.model.Product)
SpeciesContext(cbit.vcell.model.SpeciesContext)
FakeSeedSpeciesInitialConditionsParameter(org.vcell.model.rbm.FakeSeedSpeciesInitialConditionsParameter)
Reactant(cbit.vcell.model.Reactant)
BNGOutputSpec(cbit.vcell.bionetgen.BNGOutputSpec)
ExpressionException(cbit.vcell.parser.ExpressionException)
LinkedHashMap(java.util.LinkedHashMap)
FakeReactionRuleRateParameter(org.vcell.model.rbm.FakeReactionRuleRateParameter)
KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter)
RbmModelContainer(cbit.vcell.model.Model.RbmModelContainer)
Species(cbit.vcell.model.Species)
BNGSpecies(cbit.vcell.bionetgen.BNGSpecies)
HashSet(java.util.HashSet)
BNGParameter(cbit.vcell.bionetgen.BNGParameter)
ModelException(cbit.vcell.model.ModelException)
ObservableGroup(cbit.vcell.bionetgen.ObservableGroup)
RbmObservable(cbit.vcell.model.RbmObservable)
PropertyVetoException(java.beans.PropertyVetoException)
BNGReaction(cbit.vcell.bionetgen.BNGReaction)
VCUnitDefinition(cbit.vcell.units.VCUnitDefinition)
ReactionStep(cbit.vcell.model.ReactionStep)
Map(java.util.Map)
HashMap(java.util.HashMap)
LinkedHashMap(java.util.LinkedHashMap)
SpeciesPattern(org.vcell.model.rbm.SpeciesPattern)
Structure(cbit.vcell.model.Structure)
SimpleReaction(cbit.vcell.model.SimpleReaction)
ReactionRule(cbit.vcell.model.ReactionRule)
IOException(java.io.IOException)
ExpressionBindingException(cbit.vcell.parser.ExpressionBindingException)
Expression(cbit.vcell.parser.Expression)
Model(cbit.vcell.model.Model)
FakeSeedSpeciesInitialConditionsParameter(org.vcell.model.rbm.FakeSeedSpeciesInitialConditionsParameter)
Parameter(cbit.vcell.model.Parameter)
KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter)
LocalParameter(cbit.vcell.mapping.ParameterContext.LocalParameter)
BNGParameter(cbit.vcell.bionetgen.BNGParameter)
FakeReactionRuleRateParameter(org.vcell.model.rbm.FakeReactionRuleRateParameter)
MassActionKinetics(cbit.vcell.model.MassActionKinetics)
ParseException(org.vcell.model.bngl.ParseException)
BNGSpecies(cbit.vcell.bionetgen.BNGSpecies)
Example 10 with Reactant
use of cbit.vcell.model.Reactant in project vcell by virtualcell.
the class XmlReader method getReactant.
/**
* This method returns a Reactant object from a XML representation.
* Creation date: (5/4/2001 2:22:56 PM)
* @return cbit.vcell.model.Reactant
* @param param org.jdom.Element
* @exception cbit.vcell.xml.XmlParseException The exception description.
*/
private Reactant getReactant(Element param, ReactionStep reaction, Model model) throws XmlParseException {
// retrieve the key if there is one
String keystring = param.getAttributeValue(XMLTags.KeyValueAttrTag);
KeyValue key = null;
if (keystring != null && keystring.length() > 0 && this.readKeysFlag) {
key = new KeyValue(keystring);
}
String speccontref = unMangle(param.getAttributeValue(XMLTags.SpeciesContextRefAttrTag));
SpeciesContext speccont = model.getSpeciesContext(speccontref);
if (speccont == null) {
throw new XmlParseException("The reference to the SpecieContext " + speccontref + " for a SimpleReaction could not be resolved!");
}
// Retrieve Stoichiometry
int stoch = 1;
org.jdom.Attribute tempArg = param.getAttribute(XMLTags.StoichiometryAttrTag);
if (tempArg != null) {
String tempValue = tempArg.getValue();
if (tempValue.length() > 0)
stoch = Integer.parseInt(tempValue);
// param.getAttributeValue(XMLTags.StoichiometryAttrTag));
}
// return new Reactant(newkey, reaction, speccont, stoch);
return new Reactant(key, reaction, speccont, stoch);
}
Also used :
KeyValue(org.vcell.util.document.KeyValue)
Attribute(org.jdom.Attribute)
SpeciesContext(cbit.vcell.model.SpeciesContext)
Reactant(cbit.vcell.model.Reactant)
Aggregations
Reactant (cbit.vcell.model.Reactant)37
Product (cbit.vcell.model.Product)33
ReactionParticipant (cbit.vcell.model.ReactionParticipant)32
ReactionStep (cbit.vcell.model.ReactionStep)25
SpeciesContext (cbit.vcell.model.SpeciesContext)22
SimpleReaction (cbit.vcell.model.SimpleReaction)14
FluxReaction (cbit.vcell.model.FluxReaction)12
Model (cbit.vcell.model.Model)12
Structure (cbit.vcell.model.Structure)12
ArrayList (java.util.ArrayList)12
Catalyst (cbit.vcell.model.Catalyst)11
Expression (cbit.vcell.parser.Expression)10
HashMap (java.util.HashMap)9
Kinetics (cbit.vcell.model.Kinetics)8
KineticsParameter (cbit.vcell.model.Kinetics.KineticsParameter)8
ReactionRule (cbit.vcell.model.ReactionRule)8
Membrane (cbit.vcell.model.Membrane)7
Issue (org.vcell.util.Issue)7
Feature (cbit.vcell.model.Feature)6
LumpedKinetics (cbit.vcell.model.LumpedKinetics)6
