Search in sources :

Example 1 with ExperimentalPSF

use of cbit.vcell.mapping.MicroscopeMeasurement.ExperimentalPSF in project vcell by virtualcell.

the class MicroscopeMeasurementPanel method setKernel.

protected void setKernel() {
    MicroscopeMeasurement microscopeMeasurement = simulationContext.getMicroscopeMeasurement();
    if (rdbtnZprojection.isSelected()) {
        microscopeMeasurement.setConvolutionKernel(new ProjectionZKernel());
        BeanUtils.enableComponents(gaussianPsfPanel, false);
        BeanUtils.enableComponents(experimentalPsfPanel, false);
    } else if (radioButtonGaussian.isSelected()) {
        microscopeMeasurement.setConvolutionKernel(new GaussianConvolutionKernel());
        BeanUtils.enableComponents(gaussianPsfPanel, true);
        BeanUtils.enableComponents(experimentalPsfPanel, false);
    } else if (rdbtnExperimental.isSelected()) {
        String psfName = (String) pointSpreadFunctionsComboBox.getSelectedItem();
        for (DataSymbol dataSymbol : simulationContext.getDataContext().getDataSymbols()) {
            if (dataSymbol.getName().equals(psfName)) {
                microscopeMeasurement.setConvolutionKernel(new ExperimentalPSF(dataSymbol));
                break;
            }
        }
        BeanUtils.enableComponents(gaussianPsfPanel, false);
        BeanUtils.enableComponents(experimentalPsfPanel, true);
    }
}
Also used : DataSymbol(cbit.vcell.data.DataSymbol) ExperimentalPSF(cbit.vcell.mapping.MicroscopeMeasurement.ExperimentalPSF) MicroscopeMeasurement(cbit.vcell.mapping.MicroscopeMeasurement) GaussianConvolutionKernel(cbit.vcell.mapping.MicroscopeMeasurement.GaussianConvolutionKernel) ProjectionZKernel(cbit.vcell.mapping.MicroscopeMeasurement.ProjectionZKernel)

Example 2 with ExperimentalPSF

use of cbit.vcell.mapping.MicroscopeMeasurement.ExperimentalPSF 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> eventVolVarHash = new HashMap<VolVariable, EventAssignmentOrRateRuleInitParameter>();
    HashMap<Variable, RateRuleRateParameter> rateRuleRateParamHash = new HashMap<Variable, RateRuleRateParameter>();
    ArrayList<SymbolTableEntry> rateRuleVarTargets = new ArrayList<SymbolTableEntry>();
    ArrayList<SymbolTableEntry> assignmentRuleVarTargets = new ArrayList<SymbolTableEntry>();
    ArrayList<SymbolTableEntry> eventAssignTargets = 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();
        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());
                        }
                    }
                }
            }
        }
        RateRule[] rrs = simContext.getRateRules();
        if (rrs != null && rrs.length > 0) {
            for (RateRule rr : rrs) {
                SymbolTableEntry rrVar = rr.getRateRuleVar();
                if (!rateRuleVarTargets.contains(rrVar)) {
                    rateRuleVarTargets.add(rrVar);
                }
            }
        }
        AssignmentRule[] ars = simContext.getAssignmentRules();
        if (ars != null && ars.length > 0) {
            for (AssignmentRule ar : ars) {
                SymbolTableEntry arVar = ar.getAssignmentRuleVar();
                if (!assignmentRuleVarTargets.contains(arVar)) {
                    assignmentRuleVarTargets.add(arVar);
                }
            }
        }
        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])) {
            if (eventAssignTargets.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());
                }
                if (geometryClass == null) {
                    GeometryClass[] geometryClasses = simContext.getGeometryContext().getGeometry().getGeometryClasses();
                    geometryClass = geometryClasses[0];
                }
                Domain domain = null;
                if (geometryClass != null) {
                    // the volume variable will look like Compartment::g0 rather than just g0
                    domain = new Domain(geometryClass);
                }
                VolVariable volVar = new VolVariable(modelParameters[j].getName(), domain);
                varHash.addVariable(volVar);
                eventVolVarHash.put(volVar, eap);
            } else if (rateRuleVarTargets.contains(modelParameters[j])) {
                // do nothing, will do elsewhere
                ;
            } else if (assignmentRuleVarTargets.contains(modelParameters[j])) {
                // do nothing, will do elsewhere
                ;
            } else {
                Variable variable = newFunctionOrConstant(getMathSymbol(modelParameters[j], geometryClass), modelParamExpr, geometryClass);
                varHash.addVariable(variable);
            }
        }
    } 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;
                }
                String mathSymbol = getMathSymbol(parameters[i], geometryClass);
                Expression expr = getIdentifierSubstitutions(parameters[i].getExpression(), parameters[i].getUnitDefinition(), geometryClass);
                varHash.addVariable(newFunctionOrConstant(mathSymbol, expr, geometryClass));
            }
        }
    }
    // 
    // initial conditions (either function or constant) for rate rule variables that are model parameters
    // 
    // the init variables with expressions still containing variables
    Map<ModelParameter, Variable> initModelParameterHashTmp = new HashMap<>();
    // here we store the init parameter of the model parameter
    Map<EventAssignmentOrRateRuleInitParameter, ModelParameter> rateRuleInitToModelParamMapping = new HashMap<>();
    // here we store the init parameter of the model parameter
    Map<ModelParameter, EventAssignmentOrRateRuleInitParameter> modelParamTorateRuleInitMapping = new HashMap<>();
    for (ModelParameter mp : modelParameters) {
        // initial assignment for global parameter used as rate rule variable
        RateRule rr = simContext.getRateRule(mp);
        if (rr == null) {
            // we only care about global parameters that are rate rule variables
            continue;
        }
        Variable var = varHash.getVariable(mp.getName());
        if (var != null) {
            if (eventVolVarHash.containsKey(var)) {
                System.out.println("Global Parameters that are rate rule Variables should be unmapped at this point, unless they are EventAssignments too.");
            } else {
                throw new MappingException("Global Parameters that are rate rule Variables should be unmapped at this point.");
            }
        }
        Expression modelParamExpr = mp.getExpression();
        if (modelParamExpr == null) {
            continue;
        }
        GeometryClass gc = getDefaultGeometryClass(modelParamExpr);
        VCUnitDefinition paramUnit = modelUnitSystem.getInstance_TBD();
        if (mp.getUnitDefinition() != null && !mp.getUnitDefinition().equals(modelUnitSystem.getInstance_TBD())) {
            paramUnit = mp.getUnitDefinition();
        }
        // TODO: is this really needed? or could I directly use modelParamExpr in addEventAssignmentOrRateRuleInitParameter()
        Expression mpInitExpr = getIdentifierSubstitutions(modelParamExpr, paramUnit, gc);
        EventAssignmentOrRateRuleInitParameter mpInitParam;
        try {
            mpInitParam = addEventAssignmentOrRateRuleInitParameter(mp, mpInitExpr, PARAMETER_ROLE_EVENTASSIGN_OR_RATERULE_INITCONDN, paramUnit);
        } catch (PropertyVetoException e) {
            e.printStackTrace(System.out);
            throw new MappingException(e.getMessage());
        }
        rateRuleInitToModelParamMapping.put(mpInitParam, mp);
        modelParamTorateRuleInitMapping.put(mp, mpInitParam);
    }
    // 
    for (int i = 0; i < fieldMathMappingParameters.length; i++) {
        GeometryClass geometryClass = fieldMathMappingParameters[i].getGeometryClass();
        fieldMathMappingParameters[i].getExpression().bindExpression(this);
        Expression exp = getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass);
        Variable var = newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), exp, geometryClass);
        varHash.addVariable(var);
        ModelParameter mp = rateRuleInitToModelParamMapping.get(fieldMathMappingParameters[i]);
        if (mp != null) {
            initModelParameterHashTmp.put(mp, var);
        }
    }
    // 
    // initial conditions (either function or constant) for species variables
    // 
    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());
                        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 == null) {
                    String msg = initConcParm.getName() == null ? "??" : initConcParm.getName();
                    System.out.println("Unexpected NULL symbol in the initial expression of " + msg);
                } else 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());
                        initConcExpr.substituteInPlace(new Expression(spC.getName()), scsInitExpr);
                    }
                } else if (ste instanceof ModelParameter) {
                    ModelParameter mpArg = (ModelParameter) ste;
                    System.out.println(mpArg.getName());
                    if (simContext.getRateRule(mpArg) == null) {
                        // only globals that are RateRule variables need to be replaced with their _init variable
                        continue;
                    }
                    EventAssignmentOrRateRuleInitParameter mpInitParam = modelParamTorateRuleInitMapping.get(mpArg);
                    if (mpInitParam != null) {
                        // we already made it, we only need to use it
                        Expression mpArgInitExpr = new Expression(mpInitParam, getNameScope());
                        initConcExpr.substituteInPlace(new Expression(ste.getName()), mpArgInitExpr);
                    }
                } else {
                    String msg = ste.getName() == null ? "??" : ste.getName();
                    String msg2 = initConcParm.getName() == null ? "??" : initConcParm.getName();
                    System.out.println("Unexpected symbol type for " + msg + " in the initial expression of " + msg2);
                }
            }
            // 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()));
            }
        }
    }
    // 
    // functions
    // 
    enum1 = getSpeciesContextMappings();
    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 (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) {
                    // ex: eventName.delay and eventName.triggerFunction
                    String name = getMathSymbol(p, null);
                    Expression exp = getIdentifierSubstitutions(p.getExpression(), p.getUnitDefinition(), null);
                    Variable var = newFunctionOrConstant(name, exp, null);
                    varHash.addVariable(var);
                } else if (be.getParameter(BioEventParameterType.GeneralTriggerFunction) == p) {
                    // 
                    // use generated function here.
                    // 
                    String name = getMathSymbol(p, null);
                    Expression exp = getIdentifierSubstitutions(be.generateTriggerExpression(), p.getUnitDefinition(), null);
                    Variable var = newFunctionOrConstant(name, exp, null);
                    varHash.addVariable(var);
                }
            }
        }
    }
    // 
    // substitute init functions for event assignment variables
    // 
    // for (Map.Entry<VolVariable,EventAssignmentOrRateRuleInitParameter> entry : eventVolVarHash.entrySet()) {
    // EventAssignmentOrRateRuleInitParameter eap = entry.getValue();
    // 
    // String argName = eap.getName();
    // Expression modelParamExpr = eap.getExpression();
    // GeometryClass gc = getDefaultGeometryClass(modelParamExpr);
    // VCUnitDefinition paramUnit = eap.getUnitDefinition();
    // Expression mpInitExpr = new Expression(modelParamExpr);
    // String[] symbols = mpInitExpr.getSymbols();
    // if(symbols == null || symbols.length == 0) {
    // continue;
    // }
    // // TODO: this is still not working well
    // // check if 'initExpr' has other speciesContexts or rate rule global parameter variables in its expression
    // // need to replace it with 'spContext_init', modelParameter_init
    // for (String symbol : symbols) {
    // // if symbol is a speciesContext, replacing it with a reference to initial condition for that speciesContext.
    // SymbolTableEntry ste = mpInitExpr.getSymbolBinding(symbol);
    // if (ste == null) {
    // System.out.println("Unexpected NULL symbol in the initial expression of " + argName);
    // } else if (ste instanceof SpeciesContext) {
    // SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec((SpeciesContext)ste);
    // // TODO: what if initial count???
    // SpeciesContextSpecParameter spCInitParm = scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
    // // 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());
    // mpInitExpr.substituteInPlace(new Expression(ste.getName()), scsInitExpr);
    // } else if(ste instanceof ModelParameter) {
    // ModelParameter mpArg = (ModelParameter)ste;
    // System.out.println(mpArg.getName());
    // if(simContext.getRateRule(mpArg) == null) {
    // continue;		// only globals that are RateRule variables need to be replaced with their _init variable
    // }
    // Variable mpArgVar = initModelParameterHashTmp.get(mpArg);
    // if(mpArgVar != null && eventVolVarHash.get(mpArgVar) != null) {
    // EventAssignmentOrRateRuleInitParameter mpInitParam = eventVolVarHash.get(mpArgVar);
    // Expression mpArgInitExpr = new Expression(mpInitParam, getNameScope());
    // mpInitExpr.substituteInPlace(new Expression(ste.getName()), mpArgInitExpr);
    // 
    // }
    // } else {
    // String msg = ste.getName() == null ? "??" : ste.getName();
    // System.out.println("Unexpected symbol type for " + msg + " in the initial expression of " + argName);
    // }
    // }
    // varHash.removeVariable(argName);
    // Expression exp = getIdentifierSubstitutions(mpInitExpr, paramUnit, gc);
    // Variable varInit = newFunctionOrConstant(argName, exp, gc);
    // varHash.addVariable(varInit);
    // }
    // 
    // deal with rate rules
    // 
    // first, substitute the init functions for rate rule variables that are model parameters
    // we'll need this init variable (function or constant) for the ODE Equation
    // 
    // here we store the init variable with the final substitutions within their expressions
    Map<ModelParameter, Variable> initModelParameterHash = new HashMap<>();
    Map<String, SymbolTableEntry> entryMap = new HashMap<String, SymbolTableEntry>();
    simContext.getEntries(entryMap);
    for (Map.Entry<ModelParameter, Variable> entry : initModelParameterHashTmp.entrySet()) {
        ModelParameter mp = entry.getKey();
        Variable mpInitVariable = entry.getValue();
        String argName = mpInitVariable.getName();
        Expression modelParamExpr = mp.getExpression();
        GeometryClass gc = getDefaultGeometryClass(modelParamExpr);
        Expression mpInitExpr = new Expression(modelParamExpr);
        String[] symbols = mpInitExpr.getSymbols();
        if (symbols == null || symbols.length == 0) {
            // stays as it is in variable hash
            // we just move it into the initModelParameterHash
            initModelParameterHash.put(mp, mpInitVariable);
            continue;
        }
        // need to replace it with 'spContext_init', modelParameter_init
        for (String symbol : symbols) {
            // if symbol is a speciesContext, replacing it with a reference to initial condition for that speciesContext.
            SymbolTableEntry ste = mpInitExpr.getSymbolBinding(symbol);
            if (ste == null) {
                System.out.println("Unexpected NULL symbol in the initial expression of " + argName);
            } else if (ste instanceof SpeciesContext) {
                SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec((SpeciesContext) ste);
                // TODO: what if initial count???
                SpeciesContextSpecParameter spCInitParm = scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
                // 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());
                mpInitExpr.substituteInPlace(new Expression(ste.getName()), scsInitExpr);
            } else if (ste instanceof ModelParameter) {
                ModelParameter mpArg = (ModelParameter) ste;
                System.out.println(mpArg.getName());
                if (simContext.getRateRule(mpArg) == null) {
                    // only globals that are RateRule variables need to be replaced with their _init variable
                    continue;
                }
                EventAssignmentOrRateRuleInitParameter mpInitParam = modelParamTorateRuleInitMapping.get(mpArg);
                if (mpInitParam != null) {
                    // we already made it, we only need to use it
                    Expression mpArgInitExpr = new Expression(mpInitParam, getNameScope());
                    mpInitExpr.substituteInPlace(new Expression(ste.getName()), mpArgInitExpr);
                }
            } else {
                String msg = ste.getName() == null ? "??" : ste.getName();
                System.out.println("Unexpected symbol type for " + msg + " in the initial expression of " + argName);
            }
        }
        VCUnitDefinition paramUnit = modelUnitSystem.getInstance_TBD();
        if (mp.getUnitDefinition() != null && !mp.getUnitDefinition().equals(modelUnitSystem.getInstance_TBD())) {
            paramUnit = mp.getUnitDefinition();
        }
        varHash.removeVariable(mpInitVariable);
        Expression exp = getIdentifierSubstitutions(mpInitExpr, paramUnit, gc);
        mpInitVariable = newFunctionOrConstant(argName, exp, gc);
        varHash.addVariable(mpInitVariable);
        initModelParameterHash.put(mp, mpInitVariable);
    }
    // 
    // create the VolVariable for the species context used as rate rule variable
    // create the Variable (function or constant) for its rate (need it for the ODE Equation)
    // 
    enum1 = getSpeciesContextMappings();
    while (enum1.hasMoreElements()) {
        // species context used as rate rule variable
        SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
        Variable var = scm.getVariable();
        Expression exp = scm.getDependencyExpression();
        if (var == null && exp != null) {
            RateRule rr = simContext.getRateRule(scm.getSpeciesContext());
            if (rr != null && (rr.getRateRuleVar() instanceof SpeciesContext)) {
                SpeciesContext sc = scm.getSpeciesContext();
                StructureMapping sm = simContext.getGeometryContext().getStructureMapping(sc.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");
                }
                String name = getMathSymbol(scm.getSpeciesContext(), sm.getGeometryClass());
                Expression orig = rr.getRateRuleExpression();
                Expression ex = getIdentifierSubstitutions(orig, scm.getSpeciesContext().getUnitDefinition(), sm.getGeometryClass());
                GeometryClass gc = sm.getGeometryClass();
                Domain domain = null;
                if (gc != null) {
                    domain = new Domain(gc);
                }
                if (gc instanceof SurfaceClass) {
                    scm.setVariable(new MemVariable(scm.getSpeciesContext().getName(), domain));
                } else {
                    scm.setVariable(new VolVariable(scm.getSpeciesContext().getName(), domain));
                }
                Variable oldVariablre = varHash.getVariable(name);
                if (oldVariablre != null) {
                    // should always be null
                    varHash.removeVariable(name);
                }
                varHash.addVariable(scm.getVariable());
                // // create the rate parameter
                SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(sc);
                SpeciesContextSpecParameter scsInitParam = scs.getInitialConditionParameter();
                VCUnitDefinition scsInitParamUnit = scsInitParam.getUnitDefinition();
                RateRuleRateParameter rateParam = null;
                try {
                    Expression origExp = simContext.getRateRule(sc).getRateRuleExpression();
                    VCUnitDefinition rateUnit = modelUnitSystem.getInstance_TBD();
                    if (scsInitParamUnit != null && !scsInitParamUnit.equals(modelUnitSystem.getInstance_TBD())) {
                        rateUnit = scsInitParamUnit;
                    }
                    Expression rateExpr = getIdentifierSubstitutions(origExp, rateUnit, gc);
                    String argName = sc.getName() + MATH_FUNC_SUFFIX_RATERULE_RATE;
                    Variable param = newFunctionOrConstant(argName, rateExpr, gc);
                    varHash.addVariable(param);
                    rateParam = addRateRuleRateParameter(sc, rateExpr, PARAMETER_ROLE_RATERULE_RATE, rateUnit);
                } catch (PropertyVetoException e) {
                    e.printStackTrace(System.out);
                    throw new MappingException(e.getMessage());
                }
                // we generate the ODE equation elsewhere (later)
                rateRuleRateParamHash.put(scm.getVariable(), rateParam);
            }
        } else if (var != null && exp == null) {
            // could be an event variable AND a rate rule variable - in which case we need a rate parameter for the event ODE equation
            SpeciesContext sc = scm.getSpeciesContext();
            boolean isRateRuleVar = rateRuleVarTargets.contains(sc);
            boolean isEventAssignVar = eventAssignTargets.contains(sc);
            if (isRateRuleVar && isEventAssignVar) {
                // is both, so we make a rate parameter, like above
                SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(sc);
                SpeciesContextSpecParameter scsInitParam = scs.getInitialConditionParameter();
                VCUnitDefinition scsInitParamUnit = scsInitParam.getUnitDefinition();
                StructureMapping sm = simContext.getGeometryContext().getStructureMapping(sc.getStructure());
                GeometryClass gc = sm.getGeometryClass();
                RateRuleRateParameter rateParam = null;
                try {
                    Expression origExp = simContext.getRateRule(sc).getRateRuleExpression();
                    VCUnitDefinition rateUnit = modelUnitSystem.getInstance_TBD();
                    if (scsInitParamUnit != null && !scsInitParamUnit.equals(modelUnitSystem.getInstance_TBD())) {
                        rateUnit = scsInitParamUnit;
                    }
                    Expression rateExpr = getIdentifierSubstitutions(origExp, rateUnit, gc);
                    String argName = sc.getName() + MATH_FUNC_SUFFIX_RATERULE_RATE;
                    Variable param = newFunctionOrConstant(argName, rateExpr, gc);
                    varHash.addVariable(param);
                    rateParam = addRateRuleRateParameter(sc, rateExpr, PARAMETER_ROLE_RATERULE_RATE, rateUnit);
                } catch (PropertyVetoException e) {
                    e.printStackTrace(System.out);
                    throw new MappingException(e.getMessage());
                }
                // we generate the ODE equation elsewhere (later)
                rateRuleRateParamHash.put(var, rateParam);
            }
        }
    }
    // 
    for (ModelParameter mp : modelParameters) {
        // global parameter used as rate rule variable
        Variable var = varHash.getVariable(mp.getName());
        RateRule rr = simContext.getRateRule(mp);
        Expression modelParamExpr = mp.getExpression();
        if (var == null && rr != null) {
            // at this point var should be a constant
            // we're under the assumption that it's non-spatial
            GeometryClass[] geometryClasses = simContext.getGeometryContext().getGeometry().getGeometryClasses();
            GeometryClass gc = geometryClasses[0];
            // SubDomain subDomain = mathDesc.getSubDomains().nextElement();
            // GeometryClass gc = getDefaultGeometryClass(modelParamExpr);
            Domain domain = null;
            if (gc != null) {
                domain = new Domain(gc);
            }
            Variable variable;
            if (gc instanceof SurfaceClass) {
                variable = new MemVariable(mp.getName(), domain);
            } else {
                variable = new VolVariable(mp.getName(), domain);
            }
            varHash.addVariable(variable);
            RateRuleRateParameter rateParam = null;
            try {
                Expression origExp = rr.getRateRuleExpression();
                VCUnitDefinition rateUnit = modelUnitSystem.getInstance_TBD();
                if (mp.getUnitDefinition() != null && !mp.getUnitDefinition().equals(modelUnitSystem.getInstance_TBD())) {
                    rateUnit = mp.getUnitDefinition().divideBy(timeUnit);
                }
                Expression rateExpr = getIdentifierSubstitutions(origExp, rateUnit, gc);
                String argName = mp.getName() + MATH_FUNC_SUFFIX_RATERULE_RATE;
                Variable param = newFunctionOrConstant(argName, rateExpr, gc);
                varHash.addVariable(param);
                rateParam = addRateRuleRateParameter(mp, rateExpr, PARAMETER_ROLE_RATERULE_RATE, rateUnit);
            } catch (PropertyVetoException e) {
                e.printStackTrace(System.out);
                throw new MappingException(e.getMessage());
            }
            // no need to put it in the hash, we make the ODE Equation right here
            // rateRuleRateParamHash.put(variable, rateParam);
            // we know it's non-spatial
            SubDomain subDomain = mathDesc.getSubDomains().nextElement();
            Equation equation = null;
            // TODO: replace the expression with the variable  ex: "g0_protocol_init" computed above
            Expression initial = new Expression(mp.getExpression());
            // TODO: can it be null? should check and maybe try mp.getConstantValue() too ???
            Variable mpInitVariable = initModelParameterHash.get(mp);
            if (mpInitVariable != null) {
                initial = new Expression(mpInitVariable.getName());
            }
            Expression rateExpr = new Expression(0.0);
            // RateRuleRateParameter rateParam = rateRuleRateParamHash.get(variable);
            if (rateParam != null) {
                // ex: g0_rate
                rateExpr = new Expression(getMathSymbol(rateParam, gc));
            }
            // ODE Equation for rate rule variable being a global parameter
            equation = new OdeEquation(variable, initial, rateExpr);
            subDomain.addEquation(equation);
        }
    }
    // 
    // deal with assignment rules
    // 
    enum1 = getSpeciesContextMappings();
    while (enum1.hasMoreElements()) {
        // species context used as assignment rule variable
        SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
        if (scm.getVariable() == null && scm.getDependencyExpression() != null) {
            AssignmentRule ar = simContext.getAssignmentRule(scm.getSpeciesContext());
            if (ar != null && (ar.getAssignmentRuleVar() instanceof SpeciesContext)) {
                // TODO: we limit assignment rules to SpeciesContext for now
                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");
                }
                String name = getMathSymbol(scm.getSpeciesContext(), sm.getGeometryClass());
                Expression orig = ar.getAssignmentRuleExpression();
                Expression ex = getIdentifierSubstitutions(orig, scm.getSpeciesContext().getUnitDefinition(), sm.getGeometryClass());
                GeometryClass gc = sm.getGeometryClass();
                Variable dependentVariable = newFunctionOrConstant(name, ex, gc);
                dependentVariable.setDomain(new Domain(sm.getGeometryClass()));
                varHash.removeVariable(name);
                varHash.addVariable(dependentVariable);
            }
        }
    }
    for (ModelParameter mp : modelParameters) {
        // global parameter used as assignment rule variable
        Variable var = varHash.getVariable(mp.getName());
        AssignmentRule ar = simContext.getAssignmentRule(mp);
        Expression modelParamExpr = mp.getExpression();
        if (var == null && ar != null) {
            // at this point var (global parameter used as assignment rule variable) should be null
            // we're under the assumption that it's non-spatial
            GeometryClass[] geometryClasses = simContext.getGeometryContext().getGeometry().getGeometryClasses();
            GeometryClass gc = geometryClasses[0];
            SubDomain subDomain = mathDesc.getSubDomains().nextElement();
            Expression origExp = ar.getAssignmentRuleExpression();
            VCUnitDefinition rateUnit = modelUnitSystem.getInstance_TBD();
            if (mp.getUnitDefinition() != null && !mp.getUnitDefinition().equals(modelUnitSystem.getInstance_TBD())) {
                rateUnit = mp.getUnitDefinition();
            }
            Expression rateExpr = getIdentifierSubstitutions(origExp, rateUnit, gc);
            String argName = mp.getName();
            Variable param = newFunctionOrConstant(argName, rateExpr, gc);
            varHash.addVariable(param);
        }
    }
    // 
    // 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());
                        // 
                        // if it's an event assignment variable AND a rate rule variable
                        // we replace the event rate computed above (which should be zero) with the RateRuleParameter expression
                        // 
                        RateRuleRateParameter rateParam = rateRuleRateParamHash.get(variable);
                        if (rateParam != null) {
                            rate = new Expression(getMathSymbol(rateParam, null));
                        }
                        equation = new OdeEquation(variable, initial, rate);
                        subDomain.replaceEquation(equation);
                    }
                }
            } else if (scm.getVariable() instanceof VolVariable && scm.getDependencyExpression() != null) {
                // rate rule variables are like this
                RateRule rr = simContext.getRateRule(scm.getSpeciesContext());
                if (rr != null && (rr.getRateRuleVar() instanceof SpeciesContext)) {
                    // 
                    // we generate rate rule ODE equation only for species variable that are NOT event assignment variable (see right above)
                    // for global parameters variable we do it elsewhere
                    // 
                    VolVariable variable = (VolVariable) scm.getVariable();
                    Equation equation = null;
                    if (sm.getGeometryClass() == subVolume) {
                        Expression initial = new Expression(getMathSymbol(initConcParameter, null));
                        Expression rateExpr = new Expression(0.0);
                        RateRuleRateParameter rateParam = rateRuleRateParamHash.get(variable);
                        if (rateParam != null) {
                            rateExpr = new Expression(getMathSymbol(rateParam, null));
                        }
                        equation = new OdeEquation(variable, initial, rateExpr);
                        subDomain.addEquation(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 = eventVolVarHash.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 = eventVolVarHash.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) Map(java.util.Map) HashMap(java.util.HashMap) 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) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) SubDomain(cbit.vcell.math.SubDomain) PointSubDomain(cbit.vcell.math.PointSubDomain) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) PropertyVetoException(java.beans.PropertyVetoException) PdeEquation(cbit.vcell.math.PdeEquation) 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)

Example 3 with ExperimentalPSF

use of cbit.vcell.mapping.MicroscopeMeasurement.ExperimentalPSF in project vcell by virtualcell.

the class MicroscopeMeasurementPanel method refreshInterface.

protected void refreshInterface() {
    if (simulationContext == null) {
        return;
    }
    MicroscopeMeasurement microscopeMeasurement = simulationContext.getMicroscopeMeasurement();
    nameTextField.setText(microscopeMeasurement.getName());
    refreshFluorescenceSpeciesList(microscopeMeasurement);
    refreshAllSpeciesList();
    ConvolutionKernel ck = microscopeMeasurement.getConvolutionKernel();
    if (ck instanceof ProjectionZKernel) {
        rdbtnZprojection.setSelected(true);
        BeanUtils.enableComponents(gaussianPsfPanel, false);
        BeanUtils.enableComponents(experimentalPsfPanel, false);
    } else if (ck instanceof GaussianConvolutionKernel) {
        radioButtonGaussian.setSelected(true);
        BeanUtils.enableComponents(gaussianPsfPanel, true);
        BeanUtils.enableComponents(experimentalPsfPanel, false);
        sigmaXYTextField.setText(((GaussianConvolutionKernel) ck).getSigmaXY_um().infix());
        sigmaZTextField.setText(((GaussianConvolutionKernel) ck).getSigmaZ_um().infix());
    } else if (ck instanceof ExperimentalPSF) {
        rdbtnExperimental.setSelected(true);
        BeanUtils.enableComponents(gaussianPsfPanel, false);
        BeanUtils.enableComponents(experimentalPsfPanel, true);
    }
    pointSpreadFunctionsComboModel.removeAllElements();
    if (simulationContext.getDataContext() != null) {
        for (DataSymbol dataSymbol : simulationContext.getDataContext().getDataSymbols()) {
            if (dataSymbol.getDataSymbolType().equals(DataSymbolType.POINT_SPREAD_FUNCTION)) {
                pointSpreadFunctionsComboModel.addElement(dataSymbol.getName());
            }
        }
    }
}
Also used : DataSymbol(cbit.vcell.data.DataSymbol) GaussianConvolutionKernel(cbit.vcell.mapping.MicroscopeMeasurement.GaussianConvolutionKernel) ConvolutionKernel(cbit.vcell.mapping.MicroscopeMeasurement.ConvolutionKernel) ExperimentalPSF(cbit.vcell.mapping.MicroscopeMeasurement.ExperimentalPSF) MicroscopeMeasurement(cbit.vcell.mapping.MicroscopeMeasurement) GaussianConvolutionKernel(cbit.vcell.mapping.MicroscopeMeasurement.GaussianConvolutionKernel) ProjectionZKernel(cbit.vcell.mapping.MicroscopeMeasurement.ProjectionZKernel)

Aggregations

DataSymbol (cbit.vcell.data.DataSymbol)3 ExperimentalPSF (cbit.vcell.mapping.MicroscopeMeasurement.ExperimentalPSF)3 GaussianConvolutionKernel (cbit.vcell.mapping.MicroscopeMeasurement.GaussianConvolutionKernel)3 ProjectionZKernel (cbit.vcell.mapping.MicroscopeMeasurement.ProjectionZKernel)3 MicroscopeMeasurement (cbit.vcell.mapping.MicroscopeMeasurement)2 ConvolutionKernel (cbit.vcell.mapping.MicroscopeMeasurement.ConvolutionKernel)2 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 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