Examples with DistributedKinetics cbit.vcell.model.DistributedKinetics used on opensource projects

Example 1 with DistributedKinetics

use of cbit.vcell.model.DistributedKinetics in project vcell by virtualcell.

the class SBMLExporter method addReactions.

/**
 * addReactions comment.
 * @throws SbmlException
 * @throws XMLStreamException
 */
protected void addReactions() throws SbmlException, XMLStreamException {
    // Check if any reaction has electrical mapping
    boolean bCalculatePotential = false;
    StructureMapping[] structureMappings = getSelectedSimContext().getGeometryContext().getStructureMappings();
    for (int i = 0; i < structureMappings.length; i++) {
        if (structureMappings[i] instanceof MembraneMapping) {
            if (((MembraneMapping) structureMappings[i]).getCalculateVoltage()) {
                bCalculatePotential = true;
            }
        }
    }
    // If it does, VCell doesn't export it to SBML (no representation).
    if (bCalculatePotential) {
        throw new RuntimeException("This VCell model has Electrical mapping; cannot be exported to SBML at this time");
    }
    l2gMap.clear();
    ReactionSpec[] vcReactionSpecs = getSelectedSimContext().getReactionContext().getReactionSpecs();
    for (int i = 0; i < vcReactionSpecs.length; i++) {
        if (vcReactionSpecs[i].isExcluded()) {
            continue;
        }
        ReactionStep vcReactionStep = vcReactionSpecs[i].getReactionStep();
        // Create sbml reaction
        String rxnName = vcReactionStep.getName();
        org.sbml.jsbml.Reaction sbmlReaction = sbmlModel.createReaction();
        sbmlReaction.setId(org.vcell.util.TokenMangler.mangleToSName(rxnName));
        sbmlReaction.setName(rxnName);
        // If the reactionStep is a flux reaction, add the details to the annotation (structure, carrier valence, flux carrier, fluxOption, etc.)
        // If reactionStep is a simple reaction, add annotation to indicate the structure of reaction.
        // Useful when roundtripping ...
        Element sbmlImportRelatedElement = null;
        // try {
        // sbmlImportRelatedElement = getAnnotationElement(vcReactionStep);
        // } catch (XmlParseException e1) {
        // e1.printStackTrace(System.out);
        // //			throw new RuntimeException("Error ");
        // }
        // Get annotation (RDF and non-RDF) for reactionStep from SBMLAnnotationUtils
        sbmlAnnotationUtil.writeAnnotation(vcReactionStep, sbmlReaction, sbmlImportRelatedElement);
        // Now set notes,
        sbmlAnnotationUtil.writeNotes(vcReactionStep, sbmlReaction);
        // Get reaction kineticLaw
        Kinetics vcRxnKinetics = vcReactionStep.getKinetics();
        org.sbml.jsbml.KineticLaw sbmlKLaw = sbmlReaction.createKineticLaw();
        try {
            // Convert expression from kinetics rate parameter into MathML and use libSBMl utilities to convert it to formula
            // (instead of directly using rate parameter's expression infix) to maintain integrity of formula :
            // for example logical and inequalities are not handled gracefully by libSBMl if expression.infix is used.
            final Expression localRateExpr;
            final Expression lumpedRateExpr;
            if (vcRxnKinetics instanceof DistributedKinetics) {
                localRateExpr = ((DistributedKinetics) vcRxnKinetics).getReactionRateParameter().getExpression();
                lumpedRateExpr = null;
            } else if (vcRxnKinetics instanceof LumpedKinetics) {
                localRateExpr = null;
                lumpedRateExpr = ((LumpedKinetics) vcRxnKinetics).getLumpedReactionRateParameter().getExpression();
            } else {
                throw new RuntimeException("unexpected Rate Law '" + vcRxnKinetics.getClass().getSimpleName() + "', not distributed or lumped type");
            }
            // if (vcRxnKinetics instanceof DistributedKinetics)
            // Expression correctedRateExpr = kineticsAdapter.getExpression();
            // Add parameters, if any, to the kineticLaw
            Kinetics.KineticsParameter[] vcKineticsParams = vcRxnKinetics.getKineticsParameters();
            // In the first pass thro' the kinetic params, store the non-numeric param names and expressions in arrays
            String[] kinParamNames = new String[vcKineticsParams.length];
            Expression[] kinParamExprs = new Expression[vcKineticsParams.length];
            for (int j = 0; j < vcKineticsParams.length; j++) {
                if (true) {
                    // Since local reaction parameters cannot be defined by a rule, such parameters (with rules) are exported as global parameters.
                    if ((vcKineticsParams[j].getRole() == Kinetics.ROLE_CurrentDensity && (!vcKineticsParams[j].getExpression().isZero())) || (vcKineticsParams[j].getRole() == Kinetics.ROLE_LumpedCurrent && (!vcKineticsParams[j].getExpression().isZero()))) {
                        throw new RuntimeException("Electric current not handled by SBML export; failed to export reaction \"" + vcReactionStep.getName() + "\" at this time");
                    }
                    if (!vcKineticsParams[j].getExpression().isNumeric()) {
                        // NON_NUMERIC KINETIC PARAM
                        // Create new name for kinetic parameter and store it in kinParamNames, store corresponding exprs in kinParamExprs
                        // Will be used later to add this param as global.
                        String newParamName = TokenMangler.mangleToSName(vcKineticsParams[j].getName() + "_" + vcReactionStep.getName());
                        kinParamNames[j] = newParamName;
                        kinParamExprs[j] = new Expression(vcKineticsParams[j].getExpression());
                    }
                }
            }
            // If so, these need to be added as global param (else the SBML doc will not be valid)
            for (int j = 0; j < vcKineticsParams.length; j++) {
                final KineticsParameter vcKParam = vcKineticsParams[j];
                if ((vcKParam.getRole() != Kinetics.ROLE_ReactionRate) && (vcKParam.getRole() != Kinetics.ROLE_LumpedReactionRate)) {
                    // if expression of kinetic param evaluates to a double, the parameter value is set
                    if ((vcKParam.getRole() == Kinetics.ROLE_CurrentDensity && (!vcKParam.getExpression().isZero())) || (vcKParam.getRole() == Kinetics.ROLE_LumpedCurrent && (!vcKParam.getExpression().isZero()))) {
                        throw new RuntimeException("Electric current not handled by SBML export; failed to export reaction \"" + vcReactionStep.getName() + "\" at this time");
                    }
                    if (vcKParam.getExpression().isNumeric()) {
                        // NUMERIC KINETIC PARAM
                        // check if it is used in other parameters that have expressions,
                        boolean bAddedParam = false;
                        String origParamName = vcKParam.getName();
                        String newParamName = TokenMangler.mangleToSName(origParamName + "_" + vcReactionStep.getName());
                        VCUnitDefinition vcUnit = vcKParam.getUnitDefinition();
                        for (int k = 0; k < vcKineticsParams.length; k++) {
                            if (kinParamExprs[k] != null) {
                                // The param could be in the expression for any other param
                                if (kinParamExprs[k].hasSymbol(origParamName)) {
                                    // mangle its name to avoid conflict with other globals
                                    if (globalParamNamesHash.get(newParamName) == null) {
                                        globalParamNamesHash.put(newParamName, newParamName);
                                        org.sbml.jsbml.Parameter sbmlKinParam = sbmlModel.createParameter();
                                        sbmlKinParam.setId(newParamName);
                                        sbmlKinParam.setValue(vcKParam.getConstantValue());
                                        final boolean constValue = vcKParam.isConstant();
                                        sbmlKinParam.setConstant(true);
                                        // Set SBML units for sbmlParam using VC units from vcParam
                                        if (!vcUnit.isTBD()) {
                                            UnitDefinition unitDefn = getOrCreateSBMLUnit(vcUnit);
                                            sbmlKinParam.setUnits(unitDefn);
                                        }
                                        Pair<String, String> origParam = new Pair<String, String>(rxnName, origParamName);
                                        l2gMap.put(origParam, newParamName);
                                        bAddedParam = true;
                                    } else {
                                    // need to get another name for param and need to change all its refereces in the other kinParam euqations.
                                    }
                                    // update the expression to contain new name, since the globalparam has new name
                                    kinParamExprs[k].substituteInPlace(new Expression(origParamName), new Expression(newParamName));
                                }
                            }
                        }
                        // If the param hasn't been added yet, it is definitely a local param. add it to kineticLaw now.
                        if (!bAddedParam) {
                            org.sbml.jsbml.LocalParameter sbmlKinParam = sbmlKLaw.createLocalParameter();
                            sbmlKinParam.setId(origParamName);
                            sbmlKinParam.setValue(vcKParam.getConstantValue());
                            System.out.println("tis constant " + sbmlKinParam.isExplicitlySetConstant());
                            // Set SBML units for sbmlParam using VC units from vcParam
                            if (!vcUnit.isTBD()) {
                                UnitDefinition unitDefn = getOrCreateSBMLUnit(vcUnit);
                                sbmlKinParam.setUnits(unitDefn);
                            }
                        } else {
                            // hence change its occurance in rate expression if it contains that param name
                            if (localRateExpr != null && localRateExpr.hasSymbol(origParamName)) {
                                localRateExpr.substituteInPlace(new Expression(origParamName), new Expression(newParamName));
                            }
                            if (lumpedRateExpr != null && lumpedRateExpr.hasSymbol(origParamName)) {
                                lumpedRateExpr.substituteInPlace(new Expression(origParamName), new Expression(newParamName));
                            }
                        }
                    }
                }
            }
            // (using the kinParamNames and kinParamExprs above) to ensure uniqueness in the global parameter names.
            for (int j = 0; j < vcKineticsParams.length; j++) {
                if (((vcKineticsParams[j].getRole() != Kinetics.ROLE_ReactionRate) && (vcKineticsParams[j].getRole() != Kinetics.ROLE_LumpedReactionRate)) && !(vcKineticsParams[j].getExpression().isNumeric())) {
                    String oldName = vcKineticsParams[j].getName();
                    String newName = kinParamNames[j];
                    // change the name of this parameter in the rate expression
                    if (localRateExpr != null && localRateExpr.hasSymbol(oldName)) {
                        localRateExpr.substituteInPlace(new Expression(oldName), new Expression(newName));
                    }
                    if (lumpedRateExpr != null && lumpedRateExpr.hasSymbol(oldName)) {
                        lumpedRateExpr.substituteInPlace(new Expression(oldName), new Expression(newName));
                    }
                    // Change the occurence of this param in other param expressions
                    for (int k = 0; k < vcKineticsParams.length; k++) {
                        if (((vcKineticsParams[k].getRole() != Kinetics.ROLE_ReactionRate) && (vcKineticsParams[j].getRole() != Kinetics.ROLE_LumpedReactionRate)) && !(vcKineticsParams[k].getExpression().isNumeric())) {
                            if (k != j && vcKineticsParams[k].getExpression().hasSymbol(oldName)) {
                                // for all params except the current param represented by index j (whose name was changed)
                                kinParamExprs[k].substituteInPlace(new Expression(oldName), new Expression(newName));
                            }
                            if (k == j && vcKineticsParams[k].getExpression().hasSymbol(oldName)) {
                                throw new RuntimeException("A parameter cannot refer to itself in its expression");
                            }
                        }
                    }
                // end for - k
                }
            }
            // In the fifth pass thro' the kinetic params, the non-numeric params are added to the global params of the model
            for (int j = 0; j < vcKineticsParams.length; j++) {
                if (((vcKineticsParams[j].getRole() != Kinetics.ROLE_ReactionRate) && (vcKineticsParams[j].getRole() != Kinetics.ROLE_LumpedReactionRate)) && !(vcKineticsParams[j].getExpression().isNumeric())) {
                    // Now, add this param to the globalParamNamesHash and add a global parameter to the sbmlModel
                    String paramName = kinParamNames[j];
                    if (globalParamNamesHash.get(paramName) == null) {
                        globalParamNamesHash.put(paramName, paramName);
                    } else {
                    // need to get another name for param and need to change all its refereces in the other kinParam euqations.
                    }
                    Pair<String, String> origParam = new Pair<String, String>(rxnName, paramName);
                    // keeps its name but becomes a global (?)
                    l2gMap.put(origParam, paramName);
                    ASTNode paramFormulaNode = getFormulaFromExpression(kinParamExprs[j]);
                    AssignmentRule sbmlParamAssignmentRule = sbmlModel.createAssignmentRule();
                    sbmlParamAssignmentRule.setVariable(paramName);
                    sbmlParamAssignmentRule.setMath(paramFormulaNode);
                    org.sbml.jsbml.Parameter sbmlKinParam = sbmlModel.createParameter();
                    sbmlKinParam.setId(paramName);
                    if (!vcKineticsParams[j].getUnitDefinition().isTBD()) {
                        sbmlKinParam.setUnits(getOrCreateSBMLUnit(vcKineticsParams[j].getUnitDefinition()));
                    }
                    // Since the parameter is being specified by a Rule, its 'constant' field shoud be set to 'false' (default - true).
                    sbmlKinParam.setConstant(false);
                }
            }
            // end for (j) - fifth pass
            // After making all necessary adjustments to the rate expression, now set the sbmlKLaw.
            final ASTNode exprFormulaNode;
            if (lumpedRateExpr != null) {
                exprFormulaNode = getFormulaFromExpression(lumpedRateExpr);
            } else {
                if (bSpatial) {
                    exprFormulaNode = getFormulaFromExpression(localRateExpr);
                } else {
                    exprFormulaNode = getFormulaFromExpression(Expression.mult(localRateExpr, new Expression(vcReactionStep.getStructure().getName())));
                }
            }
            sbmlKLaw.setMath(exprFormulaNode);
        } catch (cbit.vcell.parser.ExpressionException e) {
            e.printStackTrace(System.out);
            throw new RuntimeException("Error getting value of parameter : " + e.getMessage());
        }
        // Add kineticLaw to sbmlReaction - not needed now, since we use sbmlRxn.createKLaw() ??
        // sbmlReaction.setKineticLaw(sbmlKLaw);
        // Add reactants, products, modifiers
        // Simple reactions have catalysts, fluxes have 'flux'
        cbit.vcell.model.ReactionParticipant[] rxnParticipants = vcReactionStep.getReactionParticipants();
        for (ReactionParticipant rxnParticpant : rxnParticipants) {
            SimpleSpeciesReference ssr = null;
            SpeciesReference sr = null;
            if (rxnParticpant instanceof cbit.vcell.model.Reactant) {
                ssr = sr = sbmlReaction.createReactant();
            } else if (rxnParticpant instanceof cbit.vcell.model.Product) {
                ssr = sr = sbmlReaction.createProduct();
            }
            if (rxnParticpant instanceof cbit.vcell.model.Catalyst) {
                ssr = sbmlReaction.createModifier();
            }
            if (ssr != null) {
                ssr.setSpecies(rxnParticpant.getSpeciesContext().getName());
            }
            if (sr != null) {
                sr.setStoichiometry(Double.parseDouble(Integer.toString(rxnParticpant.getStoichiometry())));
                String modelUniqueName = vcReactionStep.getName() + '_' + rxnParticpant.getName();
                sr.setId(TokenMangler.mangleToSName(modelUniqueName));
                // SBML-REVIEW
                sr.setConstant(true);
                // int rcode = sr.appendNotes("<
                try {
                    SBMLHelper.addNote(sr, "VCELL guess: how do we know if reaction is constant?");
                } catch (Exception e) {
                    e.printStackTrace();
                }
            }
        }
        sbmlReaction.setFast(vcReactionSpecs[i].isFast());
        // this attribute is mandatory for L3, optional for L2. So explicitly setting value.
        sbmlReaction.setReversible(true);
        if (bSpatial) {
            // set compartment for reaction if spatial
            sbmlReaction.setCompartment(vcReactionStep.getStructure().getName());
            // CORE  HAS ALT MATH true
            // set the "isLocal" attribute = true (in 'spatial' namespace) for each species
            SpatialReactionPlugin srplugin = (SpatialReactionPlugin) sbmlReaction.getPlugin(SBMLUtils.SBML_SPATIAL_NS_PREFIX);
            srplugin.setIsLocal(vcRxnKinetics instanceof DistributedKinetics);
        }
    }
}

Example 2 with DistributedKinetics

use of cbit.vcell.model.DistributedKinetics in project vcell by virtualcell.

the class StructureAnalyzer method getCorrectedRateExpression.

public Expression getCorrectedRateExpression(ReactionStep reactionStep, ReactionParticipant reactionParticipant, RateType rateType) throws Exception {
    if (reactionParticipant instanceof Catalyst) {
        throw new Exception("Catalyst " + reactionParticipant + " doesn't have a rate for this reaction");
    // return new Expression(0.0);
    }
    double stoich = reactionStep.getStoichiometry(reactionParticipant.getSpeciesContext());
    if (stoich == 0.0) {
        return new Expression(0.0);
    }
    // 
    // make distributed rate with correct stoichiometry for this participant
    // 
    VCUnitDefinition correctedReactionRateUnit = null;
    Expression distribRate = null;
    if (reactionStep.getKinetics() instanceof DistributedKinetics) {
        DistributedKinetics distributedKinetics = (DistributedKinetics) reactionStep.getKinetics();
        KineticsParameter distribReactionRateParameter = distributedKinetics.getReactionRateParameter();
        distribRate = new Expression(distribReactionRateParameter, mathMapping.getNameScope());
        correctedReactionRateUnit = distribReactionRateParameter.getUnitDefinition();
    } else if (reactionStep.getKinetics() instanceof LumpedKinetics) {
        // 
        // need to put this into concentration/time with respect to structure for reaction.
        // 
        Structure.StructureSize structureSize = reactionStep.getStructure().getStructureSize();
        LumpedKinetics lumpedKinetics = (LumpedKinetics) reactionStep.getKinetics();
        KineticsParameter lumpedReactionRateParameter = lumpedKinetics.getLumpedReactionRateParameter();
        Expression lumpedReactionRateExp = new Expression(lumpedReactionRateParameter, mathMapping.getNameScope());
        distribRate = Expression.div(lumpedReactionRateExp, new Expression(structureSize, mathMapping.getNameScope()));
        ;
        correctedReactionRateUnit = lumpedReactionRateParameter.getUnitDefinition().divideBy(structureSize.getUnitDefinition());
    }
    // correct for stoichiometry
    Expression distribRateWithStoich = distribRate;
    if (stoich != 1) {
        distribRateWithStoich = Expression.mult(new Expression(stoich), distribRateWithStoich);
    }
    // flux correction if reaction and reactionParticipant are in different compartments. (not necessarily dimensionless, use KFlux parameter).
    Expression distribRateWithStoichFlux = distribRateWithStoich;
    if (reactionStep.getStructure() != reactionParticipant.getStructure()) {
        StructureMapping reactionSM = mathMapping.getSimulationContext().getGeometryContext().getStructureMapping(reactionStep.getStructure());
        StructureMapping speciesSM = mathMapping.getSimulationContext().getGeometryContext().getStructureMapping(reactionParticipant.getStructure());
        Parameter fluxCorrectionParameter = mathMapping.getFluxCorrectionParameter(reactionSM, speciesSM);
        Expression fluxCorrection = new Expression(fluxCorrectionParameter, mathMapping.getNameScope());
        distribRateWithStoichFlux = Expression.mult(fluxCorrection, distribRateWithStoichFlux);
        correctedReactionRateUnit = correctedReactionRateUnit.multiplyBy(fluxCorrectionParameter.getUnitDefinition());
    }
    // apply unit factor for difference substance
    ModelUnitSystem unitSystem = mathMapping.getSimulationContext().getModel().getUnitSystem();
    VCUnitDefinition timeUnit = unitSystem.getTimeUnit();
    VCUnitDefinition speciesConcUnit = reactionParticipant.getSpeciesContext().getUnitDefinition();
    VCUnitDefinition speciesConcRateUnit = speciesConcUnit.divideBy(timeUnit);
    Expression unitFactor = null;
    if (rateType == RateType.ConcentrationRate) {
        unitFactor = mathMapping.getUnitFactor(speciesConcRateUnit.divideBy(correctedReactionRateUnit));
    } else if (rateType == RateType.ResolvedFluxRate) {
        unitFactor = mathMapping.getUnitFactor(speciesConcRateUnit.multiplyBy(unitSystem.getLengthUnit()).divideBy(correctedReactionRateUnit));
    }
    return Expression.mult(unitFactor, distribRateWithStoichFlux).flatten();
}

Example 3 with DistributedKinetics

use of cbit.vcell.model.DistributedKinetics in project vcell by virtualcell.

the class StructureAnalyzer method getReactionRateExpression.

public Expression getReactionRateExpression(ReactionStep reactionStep, ReactionParticipant reactionParticipant) throws Exception {
    if (reactionParticipant instanceof Catalyst) {
        throw new Exception("Catalyst " + reactionParticipant + " doesn't have a rate for this reaction");
    // return new Expression(0.0);
    }
    double stoich = reactionStep.getStoichiometry(reactionParticipant.getSpeciesContext());
    if (stoich == 0.0) {
        return new Expression(0.0);
    }
    if (reactionStep.getKinetics() instanceof DistributedKinetics) {
        DistributedKinetics distributedKinetics = (DistributedKinetics) reactionStep.getKinetics();
        if (stoich != 1) {
            Expression exp = Expression.mult(new Expression(stoich), new Expression(distributedKinetics.getReactionRateParameter(), mathMapping_4_8.getNameScope()));
            return exp;
        } else {
            Expression exp = new Expression(distributedKinetics.getReactionRateParameter(), mathMapping_4_8.getNameScope());
            return exp;
        }
    } else if (reactionStep.getKinetics() instanceof LumpedKinetics) {
        Structure.StructureSize structureSize = reactionStep.getStructure().getStructureSize();
        // 
        // need to put this into concentration/time with respect to structure for reaction.
        // 
        LumpedKinetics lumpedKinetics = (LumpedKinetics) reactionStep.getKinetics();
        Expression factor = null;
        ModelUnitSystem unitSystem = mathMapping_4_8.getSimulationContext().getModel().getUnitSystem();
        if (reactionStep.getStructure() instanceof Feature || ((reactionStep.getStructure() instanceof Membrane) && reactionStep instanceof FluxReaction)) {
            VCUnitDefinition lumpedToVolumeSubstance = unitSystem.getVolumeSubstanceUnit().divideBy(unitSystem.getLumpedReactionSubstanceUnit());
            factor = Expression.div(new Expression(lumpedToVolumeSubstance.getDimensionlessScale()), new Expression(structureSize, mathMapping_4_8.getNameScope()));
        } else if (reactionStep.getStructure() instanceof Membrane && reactionStep instanceof SimpleReaction) {
            VCUnitDefinition lumpedToVolumeSubstance = unitSystem.getMembraneSubstanceUnit().divideBy(unitSystem.getLumpedReactionSubstanceUnit());
            factor = Expression.div(new Expression(lumpedToVolumeSubstance.getDimensionlessScale()), new Expression(structureSize, mathMapping_4_8.getNameScope()));
        } else {
            throw new RuntimeException("failed to create reaction rate expression for reaction " + reactionStep.getName() + ", with kinetic type of " + reactionStep.getKinetics().getClass().getName());
        }
        if (stoich != 1) {
            Expression exp = Expression.mult(new Expression(stoich), Expression.mult(new Expression(lumpedKinetics.getLumpedReactionRateParameter(), mathMapping_4_8.getNameScope()), factor));
            return exp;
        } else {
            Expression exp = Expression.mult(new Expression(lumpedKinetics.getLumpedReactionRateParameter(), mathMapping_4_8.getNameScope()), factor);
            return exp;
        }
    } else {
        throw new RuntimeException("unexpected kinetic type " + reactionStep.getKinetics().getClass().getName());
    }
}

Example 4 with DistributedKinetics

use of cbit.vcell.model.DistributedKinetics in project vcell by virtualcell.

the class KineticsTypeTemplatePanel method getJToggleButton.

private JButton getJToggleButton() {
    if (jToggleButton == null) {
        jToggleButton = new JButton("Convert");
        jToggleButton.addActionListener(new java.awt.event.ActionListener() {

            public void actionPerformed(java.awt.event.ActionEvent e) {
                ModelUnitSystem modelUnitSystem = getKinetics().getReactionStep().getModel().getUnitSystem();
                if (getKinetics() instanceof DistributedKinetics) {
                    try {
                        reactionStep.setKinetics(LumpedKinetics.toLumpedKinetics((DistributedKinetics) getKinetics()));
                    } catch (Exception e2) {
                        e2.printStackTrace(System.out);
                        if (getKinetics().getKineticsDescription().isElectrical()) {
                            DialogUtils.showErrorDialog(KineticsTypeTemplatePanel.this, "failed to translate into General Current Kinetics [" + modelUnitSystem.getCurrentUnit().getSymbolUnicode() + "]: " + e2.getMessage(), e2);
                        } else {
                            DialogUtils.showErrorDialog(KineticsTypeTemplatePanel.this, "failed to translate into General Lumped Kinetics [" + modelUnitSystem.getLumpedReactionRateUnit().getSymbolUnicode() + "]: " + e2.getMessage(), e2);
                        }
                    }
                } else if (getKinetics() instanceof LumpedKinetics) {
                    try {
                        reactionStep.setKinetics(DistributedKinetics.toDistributedKinetics((LumpedKinetics) getKinetics()));
                    } catch (Exception e2) {
                        e2.printStackTrace(System.out);
                        if (getKinetics().getKineticsDescription().isElectrical()) {
                            DialogUtils.showErrorDialog(KineticsTypeTemplatePanel.this, "failed to translate into General Current Density Kinetics [" + modelUnitSystem.getCurrentDensityUnit().getSymbolUnicode() + "]: " + e2.getMessage(), e2);
                        } else {
                            if (getKinetics().getReactionStep().getStructure() instanceof Feature) {
                                DialogUtils.showErrorDialog(KineticsTypeTemplatePanel.this, "failed to translate into General Kinetics [" + modelUnitSystem.getVolumeReactionRateUnit().getSymbolUnicode() + "]: " + e2.getMessage(), e2);
                            } else {
                                DialogUtils.showErrorDialog(KineticsTypeTemplatePanel.this, "failed to translate into General Kinetics [" + modelUnitSystem.getMembraneReactionRateUnit().getSymbolUnicode() + "]: " + e2.getMessage(), e2);
                            }
                        }
                    }
                }
            }
        });
    }
    return jToggleButton;
}

Example 5 with DistributedKinetics

use of cbit.vcell.model.DistributedKinetics 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()));
            }
        }
    }
    // 
    // 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 ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
    }
}