Search in sources :

Example 1 with PotentialMapping

use of cbit.vcell.mapping.potential.PotentialMapping in project vcell by virtualcell.

the class DiffEquMathMapping method refreshMathDescription.

/**
 * This method was created in VisualAge.
 */
@SuppressWarnings("deprecation")
private void refreshMathDescription() throws MappingException, MatrixException, MathException, ExpressionException, ModelException {
    // All sizes must be set for new ODE models and ratios must be set for old ones.
    simContext.checkValidity();
    // 
    // temporarily place all variables in a hashtable (before binding) and discarding duplicates (check for equality)
    // 
    VariableHash varHash = new VariableHash();
    StructureMapping[] structureMappings = simContext.getGeometryContext().getStructureMappings();
    // 
    // verify that all structures are mapped to subvolumes and all subvolumes are mapped to a structure
    // 
    // Structure structures[] =
    simContext.getGeometryContext().getModel().getStructures();
    // for (int i = 0; i < structures.length; i++){
    // StructureMapping sm = simContext.getGeometryContext().getStructureMapping(structures[i]);
    // if (sm==null || (sm.getGeometryClass() == null)){
    // localIssueList.add(new Issue(structures[i], IssueCategory.StructureNotMapped,"In Application '" + simContext.getName() + "', model structure '"+structures[i].getName()+"' not mapped to a geometry subdomain",Issue.SEVERITY_WARNING));
    // }
    // }
    // SubVolume subVolumes[] = simContext.getGeometryContext().getGeometry().getGeometrySpec().getSubVolumes();
    // for (int i = 0; i < subVolumes.length; i++){
    // Structure[] mappedStructures = simContext.getGeometryContext().getStructuresFromGeometryClass(subVolumes[i]);
    // if (mappedStructures==null || mappedStructures.length==0){
    // localIssueList.add(new Issue(subVolumes[i], IssueCategory.GeometryClassNotMapped,"In Application '" + simContext.getName() + "', geometry subVolume '"+subVolumes[i].getName()+"' not mapped from a model structure",Issue.SEVERITY_WARNING));
    // }
    // }
    // deals with model parameters
    HashMap<VolVariable, EventAssignmentOrRateRuleInitParameter> eventOrRateRuleVolVarHash = new HashMap<VolVariable, EventAssignmentOrRateRuleInitParameter>();
    HashMap<VolVariable, RateRuleRateParameter> rateRuleRateParamHash = new HashMap<VolVariable, RateRuleRateParameter>();
    ArrayList<SymbolTableEntry> rateRuleVarTargets = new ArrayList<SymbolTableEntry>();
    Model model = simContext.getModel();
    ModelUnitSystem modelUnitSystem = model.getUnitSystem();
    VCUnitDefinition timeUnit = modelUnitSystem.getTimeUnit();
    ModelParameter[] modelParameters = model.getModelParameters();
    if (simContext.getGeometry().getDimension() == 0) {
        // 
        // global parameters from model (that presently are constants)
        // 
        BioEvent[] bioEvents = simContext.getBioEvents();
        ArrayList<SymbolTableEntry> eventAssignTargets = new ArrayList<SymbolTableEntry>();
        if (bioEvents != null && bioEvents.length > 0) {
            for (BioEvent be : bioEvents) {
                ArrayList<EventAssignment> eventAssignments = be.getEventAssignments();
                if (eventAssignments != null) {
                    for (EventAssignment ea : eventAssignments) {
                        if (!eventAssignTargets.contains(ea.getTarget())) {
                            eventAssignTargets.add(ea.getTarget());
                        }
                    }
                }
            }
        }
        /**
         * @author anu : RATE RULES
         */
        RateRule[] rateRules = simContext.getRateRules();
        if (rateRules != null && rateRules.length > 0) {
            for (RateRule rr : rateRules) {
                SymbolTableEntry rateRuleVar = rr.getRateRuleVar();
                if (!rateRuleVarTargets.contains(rateRuleVar)) {
                    rateRuleVarTargets.add(rateRuleVar);
                }
            }
        }
        for (int j = 0; j < modelParameters.length; j++) {
            Expression modelParamExpr = modelParameters[j].getExpression();
            GeometryClass geometryClass = getDefaultGeometryClass(modelParamExpr);
            VCUnitDefinition paramUnit = modelParameters[j].getUnitDefinition();
            modelParamExpr = getIdentifierSubstitutions(modelParamExpr, paramUnit, geometryClass);
            if (eventAssignTargets.contains(modelParameters[j]) || rateRuleVarTargets.contains(modelParameters[j])) {
                EventAssignmentOrRateRuleInitParameter eap = null;
                try {
                    eap = addEventAssignmentOrRateRuleInitParameter(modelParameters[j], modelParamExpr, PARAMETER_ROLE_EVENTASSIGN_OR_RATERULE_INITCONDN, paramUnit);
                } catch (PropertyVetoException e) {
                    e.printStackTrace(System.out);
                    throw new MappingException(e.getMessage());
                }
                VolVariable volVar = new VolVariable(modelParameters[j].getName(), null);
                varHash.addVariable(volVar);
                eventOrRateRuleVolVarHash.put(volVar, eap);
                /**
                 * RATE RULES
                 */
                if (rateRuleVarTargets.contains(modelParameters[j])) {
                    RateRuleRateParameter rateParam = null;
                    try {
                        Expression origExp = simContext.getRateRule(modelParameters[j]).getRateRuleExpression();
                        VCUnitDefinition rateUnit = modelUnitSystem.getInstance_TBD();
                        if (paramUnit != null && !paramUnit.equals(modelUnitSystem.getInstance_TBD())) {
                            rateUnit = paramUnit.divideBy(timeUnit);
                        }
                        Expression rateExpr = getIdentifierSubstitutions(origExp, rateUnit, geometryClass);
                        rateParam = addRateRuleRateParameter(modelParameters[j], rateExpr, PARAMETER_ROLE_RATERULE_RATE, rateUnit);
                    } catch (PropertyVetoException e) {
                        e.printStackTrace(System.out);
                        throw new MappingException(e.getMessage());
                    }
                    rateRuleRateParamHash.put(volVar, rateParam);
                }
            } else {
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(modelParameters[j], geometryClass), modelParamExpr, geometryClass));
            }
        }
    } else {
        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));
        }
    }
    // 
    for (SimulationContextParameter scParameter : simContext.getSimulationContextParameters()) {
        Expression scParameterExpression = scParameter.getExpression();
        GeometryClass gc = getDefaultGeometryClass(scParameterExpression);
        varHash.addVariable(newFunctionOrConstant(getMathSymbol(scParameter, gc), getIdentifierSubstitutions(scParameter.getExpression(), scParameter.getUnitDefinition(), gc), gc));
    }
    // 
    for (DataSymbol dataSymbol : simContext.getDataContext().getDataSymbols()) {
        if (dataSymbol instanceof FieldDataSymbol) {
            FieldDataSymbol fieldDataSymbol = (FieldDataSymbol) dataSymbol;
            GeometryClass geometryClass = null;
            FieldFunctionArguments ffs = new FieldFunctionArguments(fieldDataSymbol.getExternalDataIdentifier().getName(), fieldDataSymbol.getFieldDataVarName(), new Expression(fieldDataSymbol.getFieldDataVarTime()), VariableType.getVariableTypeFromVariableTypeName(fieldDataSymbol.getFieldDataVarType()));
            Expression exp = new Expression(ffs.infix());
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(dataSymbol, geometryClass), getIdentifierSubstitutions(exp, dataSymbol.getUnitDefinition(), geometryClass), geometryClass));
        } else {
            throw new RuntimeException("In Application '" + simContext.getName() + "', dataSymbol type '" + dataSymbol.getClass().getName() + "' not yet supported for math generation");
        }
    }
    // 
    // gather only those reactionSteps that are not "excluded"
    // 
    ReactionSpec[] reactionSpecs = simContext.getReactionContext().getReactionSpecs();
    Vector<ReactionStep> rsList = new Vector<ReactionStep>();
    for (int i = 0; i < reactionSpecs.length; i++) {
        if (reactionSpecs[i].isExcluded() == false) {
            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("In Application '" + simContext.getName() + "', " + reactionSteps[i].getDisplayType() + " '" + reactionSteps[i].getName() + "' contains unresolved identifier(s): " + buffer);
        }
    }
    // 
    // create new MathDescription (based on simContext's previous MathDescription if possible)
    // 
    MathDescription oldMathDesc = simContext.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(simContext.getName() + "_generated");
    }
    // 
    // volume variables
    // 
    Enumeration<SpeciesContextMapping> enum1 = getSpeciesContextMappings();
    while (enum1.hasMoreElements()) {
        SpeciesContextMapping scm = enum1.nextElement();
        if (scm.getVariable() instanceof VolVariable) {
            if (!(mathDesc.getVariable(scm.getVariable().getName()) instanceof VolVariable)) {
                varHash.addVariable(scm.getVariable());
            }
        }
    }
    // 
    // membrane variables
    // 
    enum1 = getSpeciesContextMappings();
    while (enum1.hasMoreElements()) {
        SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
        if (scm.getVariable() instanceof MemVariable) {
            varHash.addVariable(scm.getVariable());
        }
    }
    // 
    // volume region variables
    // 
    enum1 = getSpeciesContextMappings();
    while (enum1.hasMoreElements()) {
        SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
        if (scm.getVariable() instanceof VolumeRegionVariable) {
            varHash.addVariable(scm.getVariable());
        }
    }
    // 
    // membrane region variables
    // 
    enum1 = getSpeciesContextMappings();
    while (enum1.hasMoreElements()) {
        SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
        if (scm.getVariable() instanceof MembraneRegionVariable) {
            varHash.addVariable(scm.getVariable());
        }
    }
    // 
    // add compartment and membrane subdomains
    // 
    ArrayList<CompartmentSubdomainContext> compartmentSubdomainContexts = new ArrayList<CompartmentSubdomainContext>();
    ArrayList<MembraneSubdomainContext> membraneSubdomainContexts = new ArrayList<MembraneSubdomainContext>();
    addSubdomains(model, compartmentSubdomainContexts, membraneSubdomainContexts);
    // membrane velocities set on MembraneSubdomains later.
    addSpatialProcesses(varHash, compartmentSubdomainContexts, membraneSubdomainContexts);
    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(simContext.getTemperatureKelvin()), model.getTEMPERATURE().getUnitDefinition(), null)));
    // 
    // only calculate potential if at least one MembraneMapping has CalculateVoltage == true
    // 
    boolean bCalculatePotential = false;
    for (int i = 0; i < structureMappings.length; i++) {
        if (structureMappings[i] instanceof MembraneMapping) {
            if (((MembraneMapping) structureMappings[i]).getCalculateVoltage()) {
                bCalculatePotential = true;
            }
        }
    }
    potentialMapping = new PotentialMapping(simContext, this);
    if (bCalculatePotential) {
        potentialMapping.computeMath();
        // 
        // copy functions for currents and constants for capacitances
        // 
        ElectricalDevice[] devices = potentialMapping.getElectricalDevices();
        for (int j = 0; j < devices.length; j++) {
            if (devices[j] instanceof MembraneElectricalDevice) {
                MembraneElectricalDevice membraneElectricalDevice = (MembraneElectricalDevice) devices[j];
                MembraneMapping memMapping = membraneElectricalDevice.getMembraneMapping();
                Parameter specificCapacitanceParm = memMapping.getParameterFromRole(MembraneMapping.ROLE_SpecificCapacitance);
                varHash.addVariable(new Constant(getMathSymbol(specificCapacitanceParm, memMapping.getGeometryClass()), getIdentifierSubstitutions(specificCapacitanceParm.getExpression(), specificCapacitanceParm.getUnitDefinition(), memMapping.getGeometryClass())));
                ElectricalDevice.ElectricalDeviceParameter transmembraneCurrentParm = membraneElectricalDevice.getParameterFromRole(ElectricalDevice.ROLE_TransmembraneCurrent);
                ElectricalDevice.ElectricalDeviceParameter totalCurrentParm = membraneElectricalDevice.getParameterFromRole(ElectricalDevice.ROLE_TotalCurrent);
                ElectricalDevice.ElectricalDeviceParameter capacitanceParm = membraneElectricalDevice.getParameterFromRole(ElectricalDevice.ROLE_Capacitance);
                GeometryClass geometryClass = membraneElectricalDevice.getMembraneMapping().getGeometryClass();
                if (totalCurrentParm != null && /* totalCurrentDensityParm.getExpression()!=null && */
                memMapping.getCalculateVoltage()) {
                    Expression totalCurrentDensityExp = (totalCurrentParm.getExpression() != null) ? (totalCurrentParm.getExpression()) : (new Expression(0.0));
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(totalCurrentParm, geometryClass), getIdentifierSubstitutions(totalCurrentDensityExp, totalCurrentParm.getUnitDefinition(), geometryClass), geometryClass));
                }
                if (transmembraneCurrentParm != null && transmembraneCurrentParm.getExpression() != null && memMapping.getCalculateVoltage()) {
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(transmembraneCurrentParm, geometryClass), getIdentifierSubstitutions(transmembraneCurrentParm.getExpression(), transmembraneCurrentParm.getUnitDefinition(), geometryClass), geometryClass));
                }
                if (capacitanceParm != null && capacitanceParm.getExpression() != null && memMapping.getCalculateVoltage()) {
                    StructureMappingParameter sizeParameter = membraneElectricalDevice.getMembraneMapping().getSizeParameter();
                    if (simContext.getGeometry().getDimension() == 0 && (sizeParameter.getExpression() == null || sizeParameter.getExpression().isZero())) {
                        varHash.addVariable(newFunctionOrConstant(getMathSymbol(capacitanceParm, geometryClass), getIdentifierSubstitutions(Expression.mult(memMapping.getNullSizeParameterValue(), specificCapacitanceParm.getExpression()), capacitanceParm.getUnitDefinition(), geometryClass), geometryClass));
                    } else {
                        varHash.addVariable(newFunctionOrConstant(getMathSymbol(capacitanceParm, geometryClass), getIdentifierSubstitutions(capacitanceParm.getExpression(), capacitanceParm.getUnitDefinition(), geometryClass), geometryClass));
                    }
                }
                // 
                if (membraneElectricalDevice.getDependentVoltageExpression() == null) {
                    // is Voltage Independent?
                    StructureMapping.StructureMappingParameter initialVoltageParm = memMapping.getInitialVoltageParameter();
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(initialVoltageParm, memMapping.getGeometryClass()), getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), memMapping.getGeometryClass()), memMapping.getGeometryClass()));
                } else // 
                // membrane forced potential
                // 
                {
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(memMapping.getMembrane().getMembraneVoltage(), memMapping.getGeometryClass()), getIdentifierSubstitutions(membraneElectricalDevice.getDependentVoltageExpression(), memMapping.getMembrane().getMembraneVoltage().getUnitDefinition(), memMapping.getGeometryClass()), memMapping.getGeometryClass()));
                }
            } else if (devices[j] instanceof CurrentClampElectricalDevice) {
                CurrentClampElectricalDevice currentClampDevice = (CurrentClampElectricalDevice) devices[j];
                // total current = current source (no capacitance)
                Parameter totalCurrentParm = currentClampDevice.getParameterFromRole(CurrentClampElectricalDevice.ROLE_TotalCurrent);
                Parameter currentParm = currentClampDevice.getParameterFromRole(CurrentClampElectricalDevice.ROLE_TransmembraneCurrent);
                // Parameter dependentVoltage = currentClampDevice.getCurrentClampStimulus().getVoltageParameter();
                Feature deviceElectrodeFeature = currentClampDevice.getCurrentClampStimulus().getElectrode().getFeature();
                Feature groundElectrodeFeature = simContext.getGroundElectrode().getFeature();
                Membrane membrane = model.getStructureTopology().getMembrane(deviceElectrodeFeature, groundElectrodeFeature);
                GeometryClass geometryClass = null;
                if (membrane != null) {
                    StructureMapping membraneStructureMapping = simContext.getGeometryContext().getStructureMapping(membrane);
                    geometryClass = membraneStructureMapping.getGeometryClass();
                }
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(totalCurrentParm, geometryClass), getIdentifierSubstitutions(totalCurrentParm.getExpression(), totalCurrentParm.getUnitDefinition(), geometryClass), geometryClass));
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(currentParm, geometryClass), getIdentifierSubstitutions(currentParm.getExpression(), currentParm.getUnitDefinition(), geometryClass), geometryClass));
                // varHash.addVariable(newFunctionOrConstant(getMathSymbol(dependentVoltage,null),getIdentifierSubstitutions(currentClampDevice.getDependentVoltageExpression(),dependentVoltage.getUnitDefinition(),null)));
                // 
                // add user-defined parameters
                // 
                ElectricalDevice.ElectricalDeviceParameter[] parameters = currentClampDevice.getParameters();
                for (int k = 0; k < parameters.length; k++) {
                    if (parameters[k].getExpression() != null) {
                        // guards against voltage parameters that are "variable".
                        varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameters[k], null), getIdentifierSubstitutions(parameters[k].getExpression(), parameters[k].getUnitDefinition(), geometryClass), geometryClass));
                    }
                }
            } else if (devices[j] instanceof VoltageClampElectricalDevice) {
                VoltageClampElectricalDevice voltageClampDevice = (VoltageClampElectricalDevice) devices[j];
                Feature deviceElectrodeFeature = voltageClampDevice.getVoltageClampStimulus().getElectrode().getFeature();
                Feature groundElectrodeFeature = simContext.getGroundElectrode().getFeature();
                Membrane membrane = model.getStructureTopology().getMembrane(deviceElectrodeFeature, groundElectrodeFeature);
                GeometryClass geometryClass = null;
                if (membrane != null) {
                    StructureMapping membraneStructureMapping = simContext.getGeometryContext().getStructureMapping(membrane);
                    geometryClass = membraneStructureMapping.getGeometryClass();
                }
                // total current = current source (no capacitance)
                Parameter totalCurrent = voltageClampDevice.getParameterFromRole(VoltageClampElectricalDevice.ROLE_TotalCurrent);
                Parameter totalCurrentParm = voltageClampDevice.getParameterFromRole(VoltageClampElectricalDevice.ROLE_TotalCurrent);
                Parameter voltageParm = voltageClampDevice.getParameterFromRole(VoltageClampElectricalDevice.ROLE_Voltage);
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(totalCurrent, geometryClass), getIdentifierSubstitutions(totalCurrent.getExpression(), totalCurrent.getUnitDefinition(), geometryClass), geometryClass));
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(totalCurrentParm, geometryClass), getIdentifierSubstitutions(totalCurrentParm.getExpression(), totalCurrentParm.getUnitDefinition(), geometryClass), geometryClass));
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(voltageParm, geometryClass), getIdentifierSubstitutions(voltageParm.getExpression(), voltageParm.getUnitDefinition(), geometryClass), geometryClass));
                // 
                // add user-defined parameters
                // 
                ElectricalDevice.ElectricalDeviceParameter[] parameters = voltageClampDevice.getParameters();
                for (int k = 0; k < parameters.length; k++) {
                    if (parameters[k].getRole() == ElectricalDevice.ROLE_UserDefined) {
                        varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameters[k], geometryClass), getIdentifierSubstitutions(parameters[k].getExpression(), parameters[k].getUnitDefinition(), geometryClass), geometryClass));
                    }
                }
            }
        }
    } else {
        // 
        for (int j = 0; j < structureMappings.length; j++) {
            if (structureMappings[j] instanceof MembraneMapping) {
                MembraneMapping memMapping = (MembraneMapping) structureMappings[j];
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(memMapping.getMembrane().getMembraneVoltage(), memMapping.getGeometryClass()), getIdentifierSubstitutions(memMapping.getInitialVoltageParameter().getExpression(), memMapping.getInitialVoltageParameter().getUnitDefinition(), memMapping.getGeometryClass()), memMapping.getGeometryClass()));
            }
        }
    }
    // 
    for (int j = 0; j < structureMappings.length; j++) {
        if (structureMappings[j] instanceof MembraneMapping) {
            MembraneMapping membraneMapping = (MembraneMapping) structureMappings[j];
            Membrane.MembraneVoltage membraneVoltage = membraneMapping.getMembrane().getMembraneVoltage();
            ElectricalDevice[] membraneDevices = potentialMapping.getElectricalDevices(membraneMapping.getMembrane());
            // ElectricalDevice membraneDevice = null;
            for (int i = 0; i < membraneDevices.length; i++) {
                if (membraneDevices[i].hasCapacitance() && membraneDevices[i].getDependentVoltageExpression() == null) {
                    GeometryClass geometryClass = membraneMapping.getGeometryClass();
                    if (geometryClass == null) {
                        throw new MappingException("Application '" + getSimulationContext().getName() + "'\nGeometry->StructureMapping->(" + structureMappings[j].getStructure().getTypeName() + ")'" + structureMappings[j].getStructure().getName() + "' must be mapped to geometry domain.\n(see 'Problems' tab)");
                    }
                    Domain domain = new Domain(geometryClass);
                    if (membraneMapping.getCalculateVoltage() && bCalculatePotential) {
                        if (geometryClass instanceof SurfaceClass) {
                            // 
                            if (mathDesc.getVariable(Membrane.MEMBRANE_VOLTAGE_REGION_NAME) == null) {
                                // varHash.addVariable(new MembraneRegionVariable(MembraneVoltage.MEMBRANE_VOLTAGE_REGION_NAME));
                                varHash.addVariable(new MembraneRegionVariable(getMathSymbol(membraneVoltage, geometryClass), domain));
                            }
                        } else {
                            // 
                            // spatially unresolved membrane, and must solve for potential ... make VolVariable for this compartment
                            // 
                            varHash.addVariable(new VolVariable(getMathSymbol(membraneVoltage, geometryClass), domain));
                        }
                        Parameter initialVoltageParm = membraneMapping.getInitialVoltageParameter();
                        Variable initVoltageFunction = newFunctionOrConstant(getMathSymbol(initialVoltageParm, geometryClass), getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), geometryClass), geometryClass);
                        varHash.addVariable(initVoltageFunction);
                    } else {
                        // 
                        // 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);
                    }
                }
            }
        }
    }
    // 
    for (int j = 0; j < reactionSteps.length; j++) {
        ReactionStep rs = reactionSteps[j];
        if (simContext.getReactionContext().getReactionSpec(rs).isExcluded()) {
            continue;
        }
        Kinetics.KineticsParameter[] parameters = rs.getKinetics().getKineticsParameters();
        GeometryClass geometryClass = null;
        if (rs.getStructure() != null) {
            geometryClass = simContext.getGeometryContext().getStructureMapping(rs.getStructure()).getGeometryClass();
        }
        if (parameters != null) {
            for (int i = 0; i < parameters.length; i++) {
                if (((parameters[i].getRole() == Kinetics.ROLE_CurrentDensity) || (parameters[i].getRole() == Kinetics.ROLE_LumpedCurrent)) && (parameters[i].getExpression() == null || parameters[i].getExpression().isZero())) {
                    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 = simContext.getReactionContext().getSpeciesContextSpecs();
    for (int i = 0; i < speciesContextSpecs.length; i++) {
        // add initial count if present (!= null)
        SpeciesContextSpecParameter initCountParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_InitialCount);
        SpeciesContext speciesContext = speciesContextSpecs[i].getSpeciesContext();
        if (initCountParm != null && initCountParm.getExpression() != null) {
            Expression initCountExpr = new Expression(initCountParm.getExpression());
            StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContext.getStructure());
            String[] symbols = initCountExpr.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 = initCountExpr.getSymbolBinding(symbols[j]);
                if (ste instanceof SpeciesContextSpecProxyParameter) {
                    SpeciesContextSpecProxyParameter spspp = (SpeciesContextSpecProxyParameter) ste;
                    if (spspp.getTarget() instanceof SpeciesContext) {
                        spC = (SpeciesContext) spspp.getTarget();
                        SpeciesContextSpec spcspec = simContext.getReactionContext().getSpeciesContextSpec(spC);
                        SpeciesContextSpecParameter 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);
                        initCountExpr.substituteInPlace(new Expression(spC.getName()), scsInitExpr);
                    }
                }
            }
            // now create the appropriate function for the current speciesContextSpec.
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(initCountParm, sm.getGeometryClass()), getIdentifierSubstitutions(initCountExpr, initCountParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
        }
        // add initial concentration (may be derived from initial count if necessary)
        SpeciesContextSpecParameter initConcParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
        if (initConcParm != null) {
            Expression initConcExpr = null;
            if (initConcParm.getExpression() != null) {
                initConcExpr = new Expression(initConcParm.getExpression());
            } else if (initCountParm != null && initCountParm.getExpression() != null) {
                Expression structureSizeExpr = new Expression(speciesContext.getStructure().getStructureSize(), getNameScope());
                VCUnitDefinition concUnit = initConcParm.getUnitDefinition();
                VCUnitDefinition countDensityUnit = initCountParm.getUnitDefinition().divideBy(speciesContext.getStructure().getStructureSize().getUnitDefinition());
                Expression unitFactor = getUnitFactor(concUnit.divideBy(countDensityUnit));
                initConcExpr = Expression.mult(Expression.div(new Expression(initCountParm, getNameScope()), structureSizeExpr), unitFactor);
            }
            StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContext.getStructure());
            String[] symbols = initConcExpr.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 = initConcExpr.getSymbolBinding(symbols[j]);
                if (ste instanceof SpeciesContextSpecProxyParameter) {
                    SpeciesContextSpecProxyParameter spspp = (SpeciesContextSpecProxyParameter) ste;
                    if (spspp.getTarget() instanceof SpeciesContext) {
                        spC = (SpeciesContext) spspp.getTarget();
                        SpeciesContextSpec spcspec = simContext.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);
                        initConcExpr.substituteInPlace(new Expression(spC.getName()), scsInitExpr);
                    }
                }
            }
            // now create the appropriate function for the current speciesContextSpec.
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(initConcParm, sm.getGeometryClass()), getIdentifierSubstitutions(initConcExpr, initConcParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
        }
    }
    // 
    for (int i = 0; i < speciesContextSpecs.length; i++) {
        SpeciesContextMapping scm = getSpeciesContextMapping(speciesContextSpecs[i].getSpeciesContext());
        SpeciesContextSpec.SpeciesContextSpecParameter diffParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_DiffusionRate);
        if (diffParm != null && (scm.isPDERequired())) {
            StructureMapping sm = simContext.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 = simContext.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 = simContext.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(model.getKMOLE().getName(), getIdentifierSubstitutions(model.getKMOLE().getExpression(), model.getKMOLE().getUnitDefinition(), null)));
    varHash.addVariable(new Constant(model.getN_PMOLE().getName(), getIdentifierSubstitutions(model.getN_PMOLE().getExpression(), model.getN_PMOLE().getUnitDefinition(), null)));
    varHash.addVariable(new Constant(model.getKMILLIVOLTS().getName(), getIdentifierSubstitutions(model.getKMILLIVOLTS().getExpression(), model.getKMILLIVOLTS().getUnitDefinition(), null)));
    varHash.addVariable(new Constant(model.getK_GHK().getName(), getIdentifierSubstitutions(model.getK_GHK().getExpression(), model.getK_GHK().getUnitDefinition(), null)));
    // 
    for (int i = 0; i < structureMappings.length; i++) {
        StructureMapping sm = structureMappings[i];
        if (simContext.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()));
            }
        }
        StructureMappingParameter sizeParm = sm.getSizeParameter();
        if (sm.getGeometryClass() != null && sizeParm != null) {
            if (simContext.getGeometry().getDimension() == 0) {
                if (sizeParm.getExpression() != null) {
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(sizeParm, sm.getGeometryClass()), getIdentifierSubstitutions(sizeParm.getExpression(), sizeParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
                }
            } else {
                String compartmentName = sm.getGeometryClass().getName();
                VCUnitDefinition sizeUnit = sm.getSizeParameter().getUnitDefinition();
                String sizeFunctionName = null;
                if (sm instanceof MembraneMapping) {
                    MembraneMapping mm = (MembraneMapping) sm;
                    if (mm.getGeometryClass() instanceof SurfaceClass) {
                        sizeFunctionName = MathFunctionDefinitions.Function_regionArea_current.getFunctionName();
                    } else if (mm.getGeometryClass() instanceof SubVolume) {
                        sizeFunctionName = MathFunctionDefinitions.Function_regionVolume_current.getFunctionName();
                    }
                } else if (sm instanceof FeatureMapping) {
                    sizeFunctionName = MathFunctionDefinitions.Function_regionVolume_current.getFunctionName();
                } else {
                    throw new RuntimeException("structure mapping " + sm.getClass().getName() + " not yet supported");
                }
                Expression totalVolumeCorrection = sm.getStructureSizeCorrection(simContext, this);
                Expression sizeFunctionExpression = Expression.function(sizeFunctionName, new Expression[] { new Expression("'" + compartmentName + "'") });
                // sizeFunctionExpression.bindExpression(mathDesc);
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(sizeParm, sm.getGeometryClass()), getIdentifierSubstitutions(Expression.mult(totalVolumeCorrection, sizeFunctionExpression), sizeUnit, sm.getGeometryClass()), sm.getGeometryClass()));
            }
        }
    }
    // 
    for (SymbolTableEntry rrvSTE : rateRuleVarTargets) {
        if (rrvSTE instanceof SpeciesContext) {
            SpeciesContext rateRuleSpContext = (SpeciesContext) rrvSTE;
            // check if speciesContext has already been added as a vol/membrane var
            SpeciesContextMapping scm = getSpeciesContextMapping(rateRuleSpContext);
            if (scm.getVariable() instanceof VolVariable || scm.getVariable() instanceof MemVariable) {
                throw new RuntimeException("SpeciesContext '" + rrvSTE.getName() + "' has an equation and is also a rate rule variable, which is not allowed.");
            }
            // get the initial condition expression of the speciesContext.
            SpeciesContextSpec speciesContextSpec = simContext.getReactionContext().getSpeciesContextSpec(rateRuleSpContext);
            SpeciesContextSpecParameter scsInitParam = speciesContextSpec.getInitialConditionParameter();
            Expression scInitExpr = scsInitParam.getExpression();
            GeometryClass geometryClass = getDefaultGeometryClass(scInitExpr);
            VCUnitDefinition scsInitParamUnit = scsInitParam.getUnitDefinition();
            scInitExpr = getIdentifierSubstitutions(scInitExpr, scsInitParamUnit, geometryClass);
            EventAssignmentOrRateRuleInitParameter eap = null;
            try {
                // create an eventAssgnmentOrRateRuleInitParameter for the rate rule variable
                eap = addEventAssignmentOrRateRuleInitParameter(rateRuleSpContext, scInitExpr, PARAMETER_ROLE_EVENTASSIGN_OR_RATERULE_INITCONDN, scsInitParamUnit);
            } catch (PropertyVetoException e) {
                e.printStackTrace(System.out);
                throw new MappingException(e.getMessage());
            }
            // create a volVariable for this speciesContext (shouldn't have one already - since a speciesContext that has a rate rule should not participate in a reaction.
            VolVariable volVar = new VolVariable(rateRuleSpContext.getName(), null);
            varHash.addVariable(volVar);
            eventOrRateRuleVolVarHash.put(volVar, eap);
            // create the rate parameter
            RateRuleRateParameter rateParam = null;
            try {
                Expression origExp = simContext.getRateRule(rateRuleSpContext).getRateRuleExpression();
                VCUnitDefinition rateUnit = modelUnitSystem.getInstance_TBD();
                if (scsInitParamUnit != null && !scsInitParamUnit.equals(modelUnitSystem.getInstance_TBD())) {
                    rateUnit = scsInitParamUnit.divideBy(timeUnit);
                }
                Expression rateExpr = getIdentifierSubstitutions(origExp, rateUnit, geometryClass);
                rateParam = addRateRuleRateParameter(rateRuleSpContext, rateExpr, PARAMETER_ROLE_RATERULE_RATE, rateUnit);
            } catch (PropertyVetoException e) {
                e.printStackTrace(System.out);
                throw new MappingException(e.getMessage());
            }
            rateRuleRateParamHash.put(volVar, rateParam);
        }
    // end if (ste instanceof SC)
    }
    // 
    for (int i = 0; i < fieldMathMappingParameters.length; i++) {
        GeometryClass geometryClass = fieldMathMappingParameters[i].getGeometryClass();
        varHash.addVariable(newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass()));
    }
    // 
    // functions
    // 
    enum1 = getSpeciesContextMappings();
    RateRule[] rateRules = simContext.getRateRules();
    while (enum1.hasMoreElements()) {
        SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
        if (scm.getVariable() == null && scm.getDependencyExpression() != null) {
            // check if speciesContext has a rateRule; then the speciesContext should not be added as a constant
            if (rateRules == null) {
                if (simContext.getRateRule(scm.getSpeciesContext()) == null) {
                    StructureMapping sm = simContext.getGeometryContext().getStructureMapping(scm.getSpeciesContext().getStructure());
                    if (sm.getGeometryClass() == null) {
                        Structure s = sm.getStructure();
                        if (s != null) {
                            throw new RuntimeException("unmapped structure " + s.getName());
                        }
                        throw new RuntimeException("structure mapping with no structure or mapping");
                    }
                    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);
                }
            }
        }
    }
    BioEvent[] bioevents = simContext.getBioEvents();
    if (bioevents != null && bioevents.length > 0) {
        for (BioEvent be : bioevents) {
            // transform the bioEvent trigger/delay to math Event
            for (LocalParameter p : be.getEventParameters()) {
                if (p.getExpression() != null) {
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(p, null), getIdentifierSubstitutions(p.getExpression(), p.getUnitDefinition(), null), null));
                } else if (be.getParameter(BioEventParameterType.GeneralTriggerFunction) == p) {
                    // 
                    // use generated function here.
                    // 
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(p, null), getIdentifierSubstitutions(be.generateTriggerExpression(), p.getUnitDefinition(), null), null));
                }
            }
        }
    }
    // 
    // set Variables to MathDescription all at once with the order resolved by "VariableHash"
    // 
    mathDesc.setAllVariables(varHash.getAlphabeticallyOrderedVariables());
    // 
    if (simContext.getGeometryContext().getGeometry() != null) {
        try {
            mathDesc.setGeometry(simContext.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");
    }
    // 
    for (CompartmentSubdomainContext compartmentSubDomainContext : compartmentSubdomainContexts) {
        SubVolume subVolume = compartmentSubDomainContext.subvolume;
        CompartmentSubDomain subDomain = mathDesc.getCompartmentSubDomain(subVolume.getName());
        // 
        // assign boundary condition types
        // 
        StructureMapping[] mappedSMs = simContext.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) {
            if (simContext.getGeometry().getDimension() > 0) {
                subDomain.setBoundaryConditionXm(mappedFM.getBoundaryConditionTypeXm());
                subDomain.setBoundaryConditionXp(mappedFM.getBoundaryConditionTypeXp());
            }
            if (simContext.getGeometry().getDimension() > 1) {
                subDomain.setBoundaryConditionYm(mappedFM.getBoundaryConditionTypeYm());
                subDomain.setBoundaryConditionYp(mappedFM.getBoundaryConditionTypeYp());
            }
            if (simContext.getGeometry().getDimension() > 2) {
                subDomain.setBoundaryConditionZm(mappedFM.getBoundaryConditionTypeZm());
                subDomain.setBoundaryConditionZp(mappedFM.getBoundaryConditionTypeZp());
            }
        }
        // 
        // create equations
        // 
        VolumeStructureAnalyzer structureAnalyzer = getVolumeStructureAnalyzer(subVolume);
        Enumeration<SpeciesContextMapping> enumSCM = getSpeciesContextMappings();
        while (enumSCM.hasMoreElements()) {
            SpeciesContextMapping scm = enumSCM.nextElement();
            SpeciesContext sc = scm.getSpeciesContext();
            StructureMapping sm = simContext.getGeometryContext().getStructureMapping(sc.getStructure());
            SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(sc);
            // 
            // if an independent volume variable, then create equation for it (if mapped to this subDomain)
            // 
            final GeometryClass gc = sm.getGeometryClass();
            if (gc == null || !gc.getName().equals(subDomain.getName())) {
                continue;
            }
            SpeciesContextSpecParameter initConcParameter = scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
            if ((scm.getVariable() instanceof VolumeRegionVariable) && scm.getDependencyExpression() == null) {
                VolumeRegionVariable volumeRegionVariable = (VolumeRegionVariable) scm.getVariable();
                Expression initial = getIdentifierSubstitutions(new Expression(initConcParameter, getNameScope()), initConcParameter.getUnitDefinition(), sm.getGeometryClass());
                Expression rate = getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass());
                VolumeRegionEquation volumeRegionEquation = new VolumeRegionEquation(volumeRegionVariable, initial);
                volumeRegionEquation.setVolumeRateExpression(rate);
                subDomain.addEquation(volumeRegionEquation);
            } else if (scm.getVariable() instanceof VolVariable && scm.getDependencyExpression() == null) {
                VolVariable variable = (VolVariable) scm.getVariable();
                Equation equation = null;
                if (sm.getGeometryClass() == subVolume) {
                    if (scm.isPDERequired()) {
                        // 
                        // species context belongs to this subDomain
                        // 
                        Expression initial = getIdentifierSubstitutions(new Expression(initConcParameter, getNameScope()), initConcParameter.getUnitDefinition(), sm.getGeometryClass());
                        Expression rate = getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass());
                        SpeciesContextSpecParameter diffusionParameter = scs.getDiffusionParameter();
                        Expression diffusion = getIdentifierSubstitutions(new Expression(diffusionParameter, getNameScope()), diffusionParameter.getUnitDefinition(), sm.getGeometryClass());
                        equation = new PdeEquation(variable, initial, rate, diffusion);
                        ((PdeEquation) equation).setBoundaryXm((scs.getBoundaryXmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryXmParameter(), sm.getGeometryClass())));
                        ((PdeEquation) equation).setBoundaryXp((scs.getBoundaryXpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryXpParameter(), sm.getGeometryClass())));
                        ((PdeEquation) equation).setBoundaryYm((scs.getBoundaryYmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryYmParameter(), sm.getGeometryClass())));
                        ((PdeEquation) equation).setBoundaryYp((scs.getBoundaryYpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryYpParameter(), sm.getGeometryClass())));
                        ((PdeEquation) equation).setBoundaryZm((scs.getBoundaryZmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryZmParameter(), sm.getGeometryClass())));
                        ((PdeEquation) equation).setBoundaryZp((scs.getBoundaryZpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryZpParameter(), sm.getGeometryClass())));
                        if (simContext.getGeometry().getDimension() >= 1) {
                            Expression velXExp = null;
                            if (scs.getVelocityXParameter().getExpression() != null) {
                                velXExp = 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(subVolume).length;
                                    if (velX_quantities.length == 1 && numRegions == 1) {
                                        velXExp = new Expression(getMathSymbol(velX_quantities[0], sm.getGeometryClass()));
                                    } else {
                                        throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
                                    }
                                }
                            }
                            ((PdeEquation) equation).setVelocityX(velXExp);
                        }
                        if (simContext.getGeometry().getDimension() >= 2) {
                            Expression velYExp = null;
                            if (scs.getVelocityYParameter().getExpression() != null) {
                                velYExp = 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(subVolume).length;
                                    if (velY_quantities.length == 1 && numRegions == 1) {
                                        velYExp = new Expression(getMathSymbol(velY_quantities[0], sm.getGeometryClass()));
                                    } else {
                                        throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
                                    }
                                }
                            }
                            ((PdeEquation) equation).setVelocityY(velYExp);
                        }
                        if (simContext.getGeometry().getDimension() == 3) {
                            Expression velZExp = null;
                            if (scs.getVelocityZParameter().getExpression() != null) {
                                velZExp = 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(subVolume).length;
                                    if (velZ_quantities.length == 1 && numRegions == 1) {
                                        velZExp = new Expression(getMathSymbol(velZ_quantities[0], sm.getGeometryClass()));
                                    } else {
                                        throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
                                    }
                                }
                            }
                            ((PdeEquation) equation).setVelocityZ(velZExp);
                        }
                        subDomain.replaceEquation(equation);
                    } else {
                        // 
                        // ODE - species context belongs to this subDomain
                        // 
                        Expression initial = new Expression(getMathSymbol(initConcParameter, null));
                        Expression rate = (scm.getRate() == null) ? new Expression(0.0) : getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass());
                        equation = new OdeEquation(variable, initial, rate);
                        subDomain.replaceEquation(equation);
                    }
                }
            }
        }
        // 
        // create fast system (if neccessary)
        // 
        SpeciesContextMapping[] fastSpeciesContextMappings = structureAnalyzer.getFastSpeciesContextMappings();
        if (fastSpeciesContextMappings != null) {
            FastSystem fastSystem = new FastSystem(mathDesc);
            for (int i = 0; i < fastSpeciesContextMappings.length; i++) {
                SpeciesContextMapping scm = fastSpeciesContextMappings[i];
                if (scm.getFastInvariant() == null) {
                    // 
                    // independant-fast variable, create a fastRate object
                    // 
                    Expression rate = getIdentifierSubstitutions(scm.getFastRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), subVolume);
                    FastRate fastRate = new FastRate(rate);
                    fastSystem.addFastRate(fastRate);
                } else {
                    // 
                    // dependant-fast variable, create a fastInvariant object
                    // 
                    Expression rate = getIdentifierSubstitutions(scm.getFastInvariant(), modelUnitSystem.getVolumeConcentrationUnit(), subVolume);
                    FastInvariant fastInvariant = new FastInvariant(rate);
                    fastSystem.addFastInvariant(fastInvariant);
                }
            }
            subDomain.setFastSystem(fastSystem);
            // constructor calls the 'refresh' method which constructs depemdency matrix, dependent/independent vars and pseudoconstants, etc.
            // FastSystemAnalyzer fs_analyzer =
            new FastSystemAnalyzer(fastSystem, mathDesc);
        }
        // 
        // create ode's for voltages to be calculated on unresolved membranes mapped to this subVolume
        // 
        Structure[] localStructures = simContext.getGeometryContext().getStructuresFromGeometryClass(subVolume);
        for (int sIndex = 0; sIndex < localStructures.length; sIndex++) {
            if (localStructures[sIndex] instanceof Membrane) {
                Membrane membrane = (Membrane) localStructures[sIndex];
                MembraneMapping membraneMapping = (MembraneMapping) simContext.getGeometryContext().getStructureMapping(membrane);
                if ((membraneMapping.getGeometryClass() instanceof SubVolume) && membraneMapping.getCalculateVoltage()) {
                    MembraneElectricalDevice capacitiveDevice = potentialMapping.getCapacitiveDevice(membrane);
                    if (capacitiveDevice.getDependentVoltageExpression() == null) {
                        VolVariable vVar = (VolVariable) mathDesc.getVariable(getMathSymbol(capacitiveDevice.getVoltageSymbol(), membraneMapping.getGeometryClass()));
                        Expression initExp = new Expression(getMathSymbol(capacitiveDevice.getMembraneMapping().getInitialVoltageParameter(), membraneMapping.getGeometryClass()));
                        subDomain.addEquation(new OdeEquation(vVar, initExp, getIdentifierSubstitutions(potentialMapping.getOdeRHS(capacitiveDevice, this), membrane.getMembraneVoltage().getUnitDefinition().divideBy(timeUnit), membraneMapping.getGeometryClass())));
                    } else {
                    // 
                    // 
                    // 
                    }
                }
            }
        }
    }
    // 
    for (MembraneSubdomainContext memSubdomainContext : membraneSubdomainContexts) {
        MembraneSubDomain memSubDomain = memSubdomainContext.membraneSubdomain;
        SurfaceClass surfaceClass = memSubdomainContext.surfaceClass;
        for (SurfaceRegionObject surfaceRegionObject : memSubdomainContext.surfaceRegionObjects) {
            if (surfaceRegionObject.isQuantityCategoryEnabled(QuantityCategory.SurfaceVelocity)) {
                int dim = simContext.getGeometry().getDimension();
                if (dim != 2) {
                    throw new MappingException("Membrane Velocity only supported for 2D geometries");
                }
                if (simContext.getGeometry().getDimension() >= 1) {
                    SpatialQuantity velXQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.SurfaceVelocity, QuantityComponent.X);
                    Expression velXExp = new Expression(velXQuantity, simContext.getNameScope());
                    memSubDomain.setVelocityX(getIdentifierSubstitutions(velXExp, velXQuantity.getUnitDefinition(), surfaceClass));
                }
                if (simContext.getGeometry().getDimension() >= 2) {
                    SpatialQuantity velYQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.SurfaceVelocity, QuantityComponent.Y);
                    Expression velYExp = new Expression(velYQuantity, simContext.getNameScope());
                    memSubDomain.setVelocityY(getIdentifierSubstitutions(velYExp, velYQuantity.getUnitDefinition(), surfaceClass));
                }
                if (simContext.getGeometry().getDimension() == 3) {
                    SpatialQuantity velZQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.SurfaceVelocity, QuantityComponent.Z);
                    Expression velZExp = new Expression(velZQuantity, simContext.getNameScope());
                    // memSubDomain.setVelocityZ(getIdentifierSubstitutions(velZExp, velZQuantity.getUnitDefinition(), surfaceClass));
                    throw new MappingException("Membrane Velocity not supported for 2D problems");
                }
            }
        }
        // 
        // create equations for membrane-bound molecular species
        // 
        MembraneStructureAnalyzer membraneStructureAnalyzer = getMembraneStructureAnalyzer(surfaceClass);
        Enumeration<SpeciesContextMapping> enumSCM = getSpeciesContextMappings();
        while (enumSCM.hasMoreElements()) {
            SpeciesContextMapping scm = enumSCM.nextElement();
            SpeciesContext sc = scm.getSpeciesContext();
            SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(sc);
            StructureMapping sm = simContext.getGeometryContext().getStructureMapping(sc.getStructure());
            // 
            if ((scm.getVariable() instanceof MembraneRegionVariable) && scm.getDependencyExpression() == null) {
                MembraneRegionEquation equation = null;
                MembraneRegionVariable memRegionVar = (MembraneRegionVariable) scm.getVariable();
                if (sm.getGeometryClass() == surfaceClass) {
                    // 
                    // species context belongs to this subDomain
                    // 
                    Expression initial = new Expression(getMathSymbol(scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration), sm.getGeometryClass()));
                    Expression rate = getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass());
                    equation = new MembraneRegionEquation(memRegionVar, initial);
                    equation.setMembraneRateExpression(rate);
                    // equation.setUniformRateExpression(newUniformRateExpression);
                    memSubDomain.replaceEquation(equation);
                }
            } else if ((scm.getVariable() instanceof MemVariable) && scm.getDependencyExpression() == null) {
                // 
                if (sm.getGeometryClass() == surfaceClass) {
                    Equation equation = null;
                    MemVariable variable = (MemVariable) scm.getVariable();
                    if (scm.isPDERequired()) {
                        // 
                        // PDE
                        // 
                        // 
                        // species context belongs to this subDomain
                        // 
                        Expression initial = new Expression(getMathSymbol(scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration), sm.getGeometryClass()));
                        Expression rate = getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass());
                        Expression diffusion = new Expression(getMathSymbol(scs.getDiffusionParameter(), sm.getGeometryClass()));
                        equation = new PdeEquation(variable, initial, rate, diffusion);
                        ((PdeEquation) equation).setBoundaryXm((scs.getBoundaryXmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryXmParameter(), sm.getGeometryClass())));
                        ((PdeEquation) equation).setBoundaryXp((scs.getBoundaryXpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryXpParameter(), sm.getGeometryClass())));
                        ((PdeEquation) equation).setBoundaryYm((scs.getBoundaryYmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryYmParameter(), sm.getGeometryClass())));
                        ((PdeEquation) equation).setBoundaryYp((scs.getBoundaryYpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryYpParameter(), sm.getGeometryClass())));
                        ((PdeEquation) equation).setBoundaryZm((scs.getBoundaryZmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryZmParameter(), sm.getGeometryClass())));
                        ((PdeEquation) equation).setBoundaryZp((scs.getBoundaryZpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryZpParameter(), sm.getGeometryClass())));
                        memSubDomain.replaceEquation(equation);
                    } else {
                        // 
                        // ODE
                        // 
                        // 
                        // species context belongs to this subDomain
                        // 
                        Expression initial = new Expression(getMathSymbol(scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration), null));
                        Expression rate = getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass());
                        equation = new OdeEquation(variable, initial, rate);
                        memSubDomain.replaceEquation(equation);
                    }
                }
            }
        }
        Enumeration<SpeciesContextMapping> enum_scm = getSpeciesContextMappings();
        while (enum_scm.hasMoreElements()) {
            SpeciesContextMapping scm = enum_scm.nextElement();
            if (scm.isPDERequired() || scm.getVariable() instanceof VolumeRegionVariable) {
                // Species species = scm.getSpeciesContext().getSpecies();
                Variable var = scm.getVariable();
                final Domain dm = var.getDomain();
                if (dm != null) {
                    final String domainName = dm.getName();
                    if (sameName(domainName, memSubDomain.getInsideCompartment()) || sameName(domainName, memSubDomain.getOutsideCompartment())) {
                        JumpCondition jc = memSubDomain.getJumpCondition(var);
                        if (jc == null) {
                            // System.out.println("MathMapping.refreshMathDescription(), adding jump condition for diffusing variable "+var.getName()+" on membrane "+membraneStructureAnalyzer.getMembrane().getName());
                            if (var instanceof VolVariable) {
                                jc = new JumpCondition((VolVariable) var);
                            } else if (var instanceof VolumeRegionVariable) {
                                jc = new JumpCondition((VolumeRegionVariable) var);
                            } else {
                                throw new RuntimeException("unexpected Variable type " + var.getClass().getName());
                            }
                            memSubDomain.addJumpCondition(jc);
                        }
                    }
                }
            }
        }
        // 
        // set jump conditions for any volume variables or volume region variables that have explicitly defined fluxes
        // 
        ResolvedFlux[] resolvedFluxes = membraneStructureAnalyzer.getResolvedFluxes();
        if (resolvedFluxes != null) {
            for (int i = 0; i < resolvedFluxes.length; i++) {
                SpeciesContext sc = resolvedFluxes[i].getSpeciesContext();
                SpeciesContextMapping scm = getSpeciesContextMapping(sc);
                StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure());
                if (scm.getVariable() instanceof VolVariable && scm.isPDERequired()) {
                    VolVariable volVar = (VolVariable) scm.getVariable();
                    JumpCondition jc = memSubDomain.getJumpCondition(volVar);
                    if (jc == null) {
                        jc = new JumpCondition(volVar);
                        memSubDomain.addJumpCondition(jc);
                    }
                    Expression flux = getIdentifierSubstitutions(resolvedFluxes[i].getFluxExpression(), resolvedFluxes[i].getUnitDefinition(), membraneStructureAnalyzer.getSurfaceClass());
                    if (memSubDomain.getInsideCompartment().getName().equals(sm.getGeometryClass().getName())) {
                        jc.setInFlux(flux);
                    } else if (memSubDomain.getOutsideCompartment().getName().equals(sm.getGeometryClass().getName())) {
                        jc.setOutFlux(flux);
                    } else {
                        throw new RuntimeException("Application  " + simContext.getName() + " : " + scm.getSpeciesContext().getName() + " has spatially resolved flux at membrane " + scm.getSpeciesContext().getStructure().getName() + " with a non-local flux species " + scm.getSpeciesContext().getName());
                    }
                } else if (scm.getVariable() instanceof VolumeRegionVariable) {
                    VolumeRegionVariable volRegionVar = (VolumeRegionVariable) scm.getVariable();
                    JumpCondition jc = memSubDomain.getJumpCondition(volRegionVar);
                    if (jc == null) {
                        jc = new JumpCondition(volRegionVar);
                        memSubDomain.addJumpCondition(jc);
                    }
                    Expression flux = getIdentifierSubstitutions(resolvedFluxes[i].getFluxExpression(), resolvedFluxes[i].getUnitDefinition(), membraneStructureAnalyzer.getSurfaceClass());
                    if (memSubDomain.getInsideCompartment().getName().equals(sm.getGeometryClass().getName())) {
                        jc.setInFlux(flux);
                    } else if (memSubDomain.getOutsideCompartment().getName().equals(sm.getGeometryClass().getName())) {
                        jc.setOutFlux(flux);
                    } else {
                        throw new RuntimeException("Application  " + simContext.getName() + " : " + scm.getSpeciesContext().getName() + " has spatially resolved flux at membrane " + scm.getSpeciesContext().getStructure().getName() + " with a non-local flux species " + scm.getSpeciesContext().getName());
                    }
                } else {
                    throw new MappingException("Application  " + simContext.getName() + " : " + scm.getSpeciesContext().getName() + " has spatially resolved flux at membrane " + scm.getSpeciesContext().getStructure().getName() + ", but doesn't diffuse in compartment " + scm.getSpeciesContext().getStructure().getName());
                }
            }
        }
        // 
        // create fast system (if neccessary)
        // 
        SpeciesContextMapping[] fastSpeciesContextMappings = membraneStructureAnalyzer.getFastSpeciesContextMappings();
        if (fastSpeciesContextMappings != null) {
            FastSystem fastSystem = new FastSystem(mathDesc);
            for (int i = 0; i < fastSpeciesContextMappings.length; i++) {
                SpeciesContextMapping scm = fastSpeciesContextMappings[i];
                if (scm.getFastInvariant() == null) {
                    // 
                    // independant-fast variable, create a fastRate object
                    // 
                    VCUnitDefinition rateUnit = scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit);
                    FastRate fastRate = new FastRate(getIdentifierSubstitutions(scm.getFastRate(), rateUnit, surfaceClass));
                    fastSystem.addFastRate(fastRate);
                } else {
                    // 
                    // dependant-fast variable, create a fastInvariant object
                    // 
                    VCUnitDefinition invariantUnit = scm.getSpeciesContext().getUnitDefinition();
                    FastInvariant fastInvariant = new FastInvariant(getIdentifierSubstitutions(scm.getFastInvariant(), invariantUnit, surfaceClass));
                    fastSystem.addFastInvariant(fastInvariant);
                }
            }
            memSubDomain.setFastSystem(fastSystem);
            // constructor calls the 'refresh' method which constructs depemdency matrix, dependent/independent vars and pseudoconstants, etc.
            // FastSystemAnalyzer fs_analyzer =
            new FastSystemAnalyzer(fastSystem, mathDesc);
        }
        // 
        // create Membrane-region equations for potential of this resolved membrane
        // 
        Structure[] resolvedSurfaceStructures = membraneStructureAnalyzer.getStructures();
        for (int m = 0; m < resolvedSurfaceStructures.length; m++) {
            if (resolvedSurfaceStructures[m] instanceof Membrane) {
                Membrane membrane = (Membrane) resolvedSurfaceStructures[m];
                MembraneMapping membraneMapping = (MembraneMapping) simContext.getGeometryContext().getStructureMapping(membrane);
                if (membraneMapping.getCalculateVoltage()) {
                    ElectricalDevice[] membraneDevices = potentialMapping.getElectricalDevices(membrane);
                    int numCapacitiveDevices = 0;
                    MembraneElectricalDevice capacitiveDevice = null;
                    for (int i = 0; i < membraneDevices.length; i++) {
                        if (membraneDevices[i] instanceof MembraneElectricalDevice) {
                            numCapacitiveDevices++;
                            capacitiveDevice = (MembraneElectricalDevice) membraneDevices[i];
                        }
                    }
                    if (numCapacitiveDevices != 1) {
                        throw new MappingException("expecting 1 capacitive electrical device on graph edge for membrane " + membrane.getName() + ", found '" + numCapacitiveDevices + "'");
                    }
                    if (mathDesc.getVariable(getMathSymbol(capacitiveDevice.getVoltageSymbol(), membraneMapping.getGeometryClass())) instanceof MembraneRegionVariable) {
                        MembraneRegionVariable vVar = (MembraneRegionVariable) mathDesc.getVariable(getMathSymbol(capacitiveDevice.getVoltageSymbol(), membraneMapping.getGeometryClass()));
                        Parameter initialVoltageParm = capacitiveDevice.getMembraneMapping().getInitialVoltageParameter();
                        Expression initExp = getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), capacitiveDevice.getMembraneMapping().getGeometryClass());
                        MembraneRegionEquation vEquation = new MembraneRegionEquation(vVar, initExp);
                        vEquation.setMembraneRateExpression(getIdentifierSubstitutions(potentialMapping.getOdeRHS(capacitiveDevice, this), membrane.getMembraneVoltage().getUnitDefinition().divideBy(timeUnit), capacitiveDevice.getMembraneMapping().getGeometryClass()));
                        memSubDomain.addEquation(vEquation);
                    }
                }
            }
        }
    }
    // create equations for event assignment or rate rule targets that are model params/species, etc.
    Set<VolVariable> hashKeySet = eventOrRateRuleVolVarHash.keySet();
    Iterator<VolVariable> volVarsIter = hashKeySet.iterator();
    // working under the assumption that we are dealing with non-spatial math, hence only one compartment domain!
    SubDomain subDomain = mathDesc.getSubDomains().nextElement();
    while (volVarsIter.hasNext()) {
        VolVariable volVar = volVarsIter.next();
        EventAssignmentOrRateRuleInitParameter initParam = eventOrRateRuleVolVarHash.get(volVar);
        // check event initial condition, it shouldn't contain vars, we have to do it here, coz we want to substitute functions...etc.
        Expression eapExp = MathUtilities.substituteFunctions(initParam.getExpression(), mathDesc);
        if (eapExp.getSymbols() != null) {
            for (String symbol : eapExp.getSymbols()) {
                SymbolTableEntry ste = eapExp.getSymbolBinding(symbol);
                if (ste instanceof VolVariable || ste instanceof MemVariable) {
                    throw new MathException("Variables are not allowed in Event assignment initial condition.\nEvent assignment target: " + volVar.getName() + " has variable (" + symbol + ") in its expression.");
                }
            }
        }
        Expression rateExpr = new Expression(0.0);
        RateRuleRateParameter rateParam = rateRuleRateParamHash.get(volVar);
        if (rateParam != null) {
            // this is a rate rule, get its expression.
            rateExpr = new Expression(getMathSymbol(rateParam, null));
        }
        Equation equation = new OdeEquation(volVar, new Expression(getMathSymbol(initParam, null)), rateExpr);
        subDomain.addEquation(equation);
    }
    // events - add events to math desc for event assignments that have parameters as target variables
    if (bioevents != null && bioevents.length > 0) {
        for (BioEvent be : bioevents) {
            // transform the bioEvent trigger/delay to math Event
            LocalParameter genTriggerParam = be.getParameter(BioEventParameterType.GeneralTriggerFunction);
            Expression mathTriggerExpr = getIdentifierSubstitutions(new Expression(genTriggerParam, be.getNameScope()), modelUnitSystem.getInstance_DIMENSIONLESS(), null);
            Delay mathDelay = null;
            LocalParameter delayParam = be.getParameter(BioEventParameterType.TriggerDelay);
            if (delayParam != null && delayParam.getExpression() != null && !delayParam.getExpression().compareEqual(new Expression(0.0))) {
                boolean bUseValsFromTriggerTime = be.getUseValuesFromTriggerTime();
                Expression mathDelayExpr = getIdentifierSubstitutions(new Expression(delayParam, be.getNameScope()), timeUnit, null);
                mathDelay = new Delay(bUseValsFromTriggerTime, mathDelayExpr);
            }
            // now deal with (bio)event Assignment translation to math EventAssignment
            ArrayList<EventAssignment> eventAssignments = be.getEventAssignments();
            ArrayList<Event.EventAssignment> mathEventAssignmentsList = new ArrayList<Event.EventAssignment>();
            if (eventAssignments != null) {
                for (EventAssignment ea : eventAssignments) {
                    SymbolTableEntry ste = simContext.getEntry(ea.getTarget().getName());
                    if (ste instanceof StructureSize) {
                        throw new RuntimeException("Event Assignment Variable for compartment size is not supported yet");
                    }
                    VCUnitDefinition eventAssignVarUnit = ste.getUnitDefinition();
                    Variable variable = varHash.getVariable(ste.getName());
                    Event.EventAssignment mathEA = new Event.EventAssignment(variable, getIdentifierSubstitutions(ea.getAssignmentExpression(), eventAssignVarUnit, null));
                    mathEventAssignmentsList.add(mathEA);
                }
            }
            // use the translated trigger, delay and event assignments to create (math) event
            Event mathEvent = new Event(be.getName(), mathTriggerExpr, mathDelay, mathEventAssignmentsList);
            mathDesc.addEvent(mathEvent);
        }
    }
    if (simContext.getMicroscopeMeasurement() != null && simContext.getMicroscopeMeasurement().getFluorescentSpecies().size() > 0) {
        MicroscopeMeasurement measurement = simContext.getMicroscopeMeasurement();
        Expression volumeConcExp = new Expression(0.0);
        Expression membraneDensityExp = new Expression(0.0);
        for (SpeciesContext speciesContext : measurement.getFluorescentSpecies()) {
            GeometryClass geometryClass = simContext.getGeometryContext().getStructureMapping(speciesContext.getStructure()).getGeometryClass();
            StructureMapping structureMapping = simContext.getGeometryContext().getStructureMapping(speciesContext.getStructure());
            StructureMappingParameter unitSizeParameter = structureMapping.getUnitSizeParameter();
            Expression mappedSpeciesContextExpression = Expression.mult(unitSizeParameter.getExpression(), new Expression(getMathSymbol(speciesContext, geometryClass)));
            VCUnitDefinition mappedSpeciesContextUnit = unitSizeParameter.getUnitDefinition().multiplyBy(speciesContext.getUnitDefinition());
            if (geometryClass instanceof SubVolume) {
                // volume function
                int dimension = 3;
                VCUnitDefinition desiredConcUnits = model.getUnitSystem().getInstance("molecules").divideBy(model.getUnitSystem().getLengthUnit().raiseTo(new ucar.units_vcell.RationalNumber(dimension)));
                Expression unitFactor = getUnitFactor(desiredConcUnits.divideBy(mappedSpeciesContextUnit));
                volumeConcExp = Expression.add(volumeConcExp, Expression.mult(unitFactor, mappedSpeciesContextExpression)).flatten();
            } else if (geometryClass instanceof SurfaceClass) {
                // membrane function
                int dimension = 2;
                VCUnitDefinition desiredSurfaceDensityUnits = model.getUnitSystem().getInstance("molecules").divideBy(model.getUnitSystem().getLengthUnit().raiseTo(new ucar.units_vcell.RationalNumber(dimension)));
                Expression unitFactor = getUnitFactor(desiredSurfaceDensityUnits.divideBy(mappedSpeciesContextUnit));
                membraneDensityExp = Expression.add(membraneDensityExp, Expression.mult(unitFactor, mappedSpeciesContextExpression)).flatten();
            } else {
                throw new MathException("unsupported geometry mapping for microscopy measurement");
            }
        }
        ConvolutionKernel kernel = measurement.getConvolutionKernel();
        if (kernel instanceof ExperimentalPSF) {
            if (!membraneDensityExp.isZero()) {
                throw new MappingException("membrane variables and functions not yet supported for Z projection in Microcopy Measurements");
            }
            ExperimentalPSF psf = (ExperimentalPSF) kernel;
            DataSymbol psfDataSymbol = psf.getPSFDataSymbol();
            if (psfDataSymbol instanceof FieldDataSymbol) {
                FieldDataSymbol fieldDataSymbol = (FieldDataSymbol) psfDataSymbol;
                String fieldDataName = ((FieldDataSymbol) psfDataSymbol).getExternalDataIdentifier().getName();
                Expression psfExp = Expression.function(FieldFunctionDefinition.FUNCTION_name, new Expression("'" + fieldDataName + "'"), new Expression("'" + fieldDataSymbol.getFieldDataVarName() + "'"), new Expression(fieldDataSymbol.getFieldDataVarTime()), new Expression("'" + fieldDataSymbol.getFieldDataVarType() + "'"));
                varHash.addVariable(new Function("__PSF__", psfExp, null));
            }
            Expression convExp = Expression.function(ConvFunctionDefinition.FUNCTION_name, volumeConcExp, new Expression("__PSF__"));
            varHash.addVariable(newFunctionOrConstant(measurement.getName(), convExp, null));
        } else if (kernel instanceof GaussianConvolutionKernel) {
            GaussianConvolutionKernel gaussianConvolutionKernel = (GaussianConvolutionKernel) kernel;
            GaussianConvolutionDataGeneratorKernel mathKernel = new GaussianConvolutionDataGeneratorKernel(gaussianConvolutionKernel.getSigmaXY_um(), gaussianConvolutionKernel.getSigmaZ_um());
            ConvolutionDataGenerator dataGenerator = new ConvolutionDataGenerator(measurement.getName(), mathKernel, volumeConcExp, membraneDensityExp);
            mathDesc.getPostProcessingBlock().addDataGenerator(dataGenerator);
        } else if (kernel instanceof ProjectionZKernel) {
            if (mathDesc.getGeometry().getDimension() == 3) {
                if (!membraneDensityExp.isZero()) {
                    throw new MappingException("membrane variables and functions not yet supported for Z projection in Microcopy Measurements");
                }
                ProjectionDataGenerator dataGenerator = new ProjectionDataGenerator(measurement.getName(), null, ProjectionDataGenerator.Axis.z, ProjectionDataGenerator.Operation.sum, volumeConcExp);
                mathDesc.getPostProcessingBlock().addDataGenerator(dataGenerator);
            } else {
                throw new MappingException("Z Projection is only supported in 3D spatial applications.");
            }
        }
    }
    // 
    for (int i = 0; i < fieldMathMappingParameters.length; i++) {
        if (fieldMathMappingParameters[i] instanceof UnitFactorParameter) {
            GeometryClass geometryClass = fieldMathMappingParameters[i].getGeometryClass();
            Variable variable = newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass());
            if (mathDesc.getVariable(variable.getName()) == null) {
                mathDesc.addVariable(variable);
            }
        }
    }
    if (!mathDesc.isValid()) {
        System.out.println(mathDesc.getVCML_database());
        throw new MappingException("generated an invalid mathDescription: " + mathDesc.getWarning());
    }
// System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string begin ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
// System.out.println(mathDesc.getVCML());
// System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string end ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
}
Also used : MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) ArrayList(java.util.ArrayList) StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter) SpeciesContext(cbit.vcell.model.SpeciesContext) Feature(cbit.vcell.model.Feature) MemVariable(cbit.vcell.math.MemVariable) SubVolume(cbit.vcell.geometry.SubVolume) Vector(java.util.Vector) SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter) FastInvariant(cbit.vcell.math.FastInvariant) GaussianConvolutionKernel(cbit.vcell.mapping.MicroscopeMeasurement.GaussianConvolutionKernel) PropertyVetoException(java.beans.PropertyVetoException) FieldDataSymbol(cbit.vcell.data.FieldDataSymbol) DataSymbol(cbit.vcell.data.DataSymbol) VCUnitDefinition(cbit.vcell.units.VCUnitDefinition) FastSystem(cbit.vcell.math.FastSystem) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) ReactionStep(cbit.vcell.model.ReactionStep) MembraneRegionEquation(cbit.vcell.math.MembraneRegionEquation) SurfaceClass(cbit.vcell.geometry.SurfaceClass) VariableHash(cbit.vcell.math.VariableHash) ExperimentalPSF(cbit.vcell.mapping.MicroscopeMeasurement.ExperimentalPSF) ConvolutionDataGenerator(cbit.vcell.math.ConvolutionDataGenerator) GaussianConvolutionDataGeneratorKernel(cbit.vcell.math.ConvolutionDataGenerator.GaussianConvolutionDataGeneratorKernel) Structure(cbit.vcell.model.Structure) ModelUnitSystem(cbit.vcell.model.ModelUnitSystem) LocalParameter(cbit.vcell.mapping.ParameterContext.LocalParameter) VoltageClampElectricalDevice(cbit.vcell.mapping.potential.VoltageClampElectricalDevice) StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter) Expression(cbit.vcell.parser.Expression) Model(cbit.vcell.model.Model) StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter) Parameter(cbit.vcell.model.Parameter) SpeciesContextSpecProxyParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecProxyParameter) SimulationContextParameter(cbit.vcell.mapping.SimulationContext.SimulationContextParameter) ModelParameter(cbit.vcell.model.Model.ModelParameter) SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter) LocalParameter(cbit.vcell.mapping.ParameterContext.LocalParameter) Event(cbit.vcell.math.Event) ProjectionDataGenerator(cbit.vcell.math.ProjectionDataGenerator) SurfaceRegionObject(cbit.vcell.mapping.spatial.SurfaceRegionObject) FieldDataSymbol(cbit.vcell.data.FieldDataSymbol) MembraneRegionVariable(cbit.vcell.math.MembraneRegionVariable) HashMap(java.util.HashMap) MathDescription(cbit.vcell.math.MathDescription) MembraneElectricalDevice(cbit.vcell.mapping.potential.MembraneElectricalDevice) Delay(cbit.vcell.math.Event.Delay) PropertyVetoException(java.beans.PropertyVetoException) PdeEquation(cbit.vcell.math.PdeEquation) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) SubDomain(cbit.vcell.math.SubDomain) PointSubDomain(cbit.vcell.math.PointSubDomain) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) CurrentClampElectricalDevice(cbit.vcell.mapping.potential.CurrentClampElectricalDevice) SpatialQuantity(cbit.vcell.mapping.spatial.SpatialObject.SpatialQuantity) SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter) EventAssignment(cbit.vcell.mapping.BioEvent.EventAssignment) CurrentClampElectricalDevice(cbit.vcell.mapping.potential.CurrentClampElectricalDevice) MembraneElectricalDevice(cbit.vcell.mapping.potential.MembraneElectricalDevice) ElectricalDevice(cbit.vcell.mapping.potential.ElectricalDevice) VoltageClampElectricalDevice(cbit.vcell.mapping.potential.VoltageClampElectricalDevice) VolVariable(cbit.vcell.math.VolVariable) StructureSize(cbit.vcell.model.Structure.StructureSize) ProjectionZKernel(cbit.vcell.mapping.MicroscopeMeasurement.ProjectionZKernel) ModelParameter(cbit.vcell.model.Model.ModelParameter) OdeEquation(cbit.vcell.math.OdeEquation) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) SubDomain(cbit.vcell.math.SubDomain) PointSubDomain(cbit.vcell.math.PointSubDomain) Domain(cbit.vcell.math.Variable.Domain) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) JumpCondition(cbit.vcell.math.JumpCondition) GeometryClass(cbit.vcell.geometry.GeometryClass) VolVariable(cbit.vcell.math.VolVariable) MembraneRegionVariable(cbit.vcell.math.MembraneRegionVariable) PointVariable(cbit.vcell.math.PointVariable) VolumeRegionVariable(cbit.vcell.math.VolumeRegionVariable) MemVariable(cbit.vcell.math.MemVariable) Variable(cbit.vcell.math.Variable) SpeciesContextSpecProxyParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecProxyParameter) Constant(cbit.vcell.math.Constant) VolumeRegionVariable(cbit.vcell.math.VolumeRegionVariable) Function(cbit.vcell.math.Function) SymbolTableEntry(cbit.vcell.parser.SymbolTableEntry) Membrane(cbit.vcell.model.Membrane) VolumeRegionEquation(cbit.vcell.math.VolumeRegionEquation) PotentialMapping(cbit.vcell.mapping.potential.PotentialMapping) EventAssignment(cbit.vcell.mapping.BioEvent.EventAssignment) FieldFunctionArguments(cbit.vcell.field.FieldFunctionArguments) PdeEquation(cbit.vcell.math.PdeEquation) ComputeMembraneMetricEquation(cbit.vcell.math.ComputeMembraneMetricEquation) VolumeRegionEquation(cbit.vcell.math.VolumeRegionEquation) OdeEquation(cbit.vcell.math.OdeEquation) MembraneRegionEquation(cbit.vcell.math.MembraneRegionEquation) Equation(cbit.vcell.math.Equation) FastRate(cbit.vcell.math.FastRate) SimulationContextParameter(cbit.vcell.mapping.SimulationContext.SimulationContextParameter) ConvolutionKernel(cbit.vcell.mapping.MicroscopeMeasurement.ConvolutionKernel) GaussianConvolutionKernel(cbit.vcell.mapping.MicroscopeMeasurement.GaussianConvolutionKernel) MathException(cbit.vcell.math.MathException)

Aggregations

DataSymbol (cbit.vcell.data.DataSymbol)1 FieldDataSymbol (cbit.vcell.data.FieldDataSymbol)1 FieldFunctionArguments (cbit.vcell.field.FieldFunctionArguments)1 GeometryClass (cbit.vcell.geometry.GeometryClass)1 SubVolume (cbit.vcell.geometry.SubVolume)1 SurfaceClass (cbit.vcell.geometry.SurfaceClass)1 EventAssignment (cbit.vcell.mapping.BioEvent.EventAssignment)1 ConvolutionKernel (cbit.vcell.mapping.MicroscopeMeasurement.ConvolutionKernel)1 ExperimentalPSF (cbit.vcell.mapping.MicroscopeMeasurement.ExperimentalPSF)1 GaussianConvolutionKernel (cbit.vcell.mapping.MicroscopeMeasurement.GaussianConvolutionKernel)1 ProjectionZKernel (cbit.vcell.mapping.MicroscopeMeasurement.ProjectionZKernel)1 LocalParameter (cbit.vcell.mapping.ParameterContext.LocalParameter)1 SimulationContextParameter (cbit.vcell.mapping.SimulationContext.SimulationContextParameter)1 SpeciesContextSpecParameter (cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter)1 SpeciesContextSpecProxyParameter (cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecProxyParameter)1 StructureMappingParameter (cbit.vcell.mapping.StructureMapping.StructureMappingParameter)1 CurrentClampElectricalDevice (cbit.vcell.mapping.potential.CurrentClampElectricalDevice)1 ElectricalDevice (cbit.vcell.mapping.potential.ElectricalDevice)1 MembraneElectricalDevice (cbit.vcell.mapping.potential.MembraneElectricalDevice)1 PotentialMapping (cbit.vcell.mapping.potential.PotentialMapping)1