diff --git a/src/equations_set_constants.F90 b/src/equations_set_constants.F90
index 53a9bbac..de76a6a6 100644
--- a/src/equations_set_constants.F90
+++ b/src/equations_set_constants.F90
@@ -163,8 +163,12 @@ MODULE EquationsSetConstants
   INTEGER(INTG), PARAMETER :: EQUATIONS_SET_RATE_BASED_GROWTH_MODEL_SUBTYPE=39
   INTEGER(INTG), PARAMETER :: EQUATIONS_SET_COMPRESSIBLE_RATE_BASED_GROWTH_MODEL_SUBTYPE=40
   INTEGER(INTG), PARAMETER :: EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE=41
-  INTEGER(INTG), PARAMETER :: EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE=42
-  INTEGER(INTG), PARAMETER :: EQUATIONS_SET_MR_AND_GROWTH_LAW_IN_CELLML_SUBTYPE=43
+  INTEGER(INTG), PARAMETER :: EQUATIONS_SET_REFERENCE_STATE_ISO_EXPONENTIAL_SUBTYPE=42
+
+  INTEGER(INTG), PARAMETER :: EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE=43
+  INTEGER(INTG), PARAMETER :: EQUATIONS_SET_MR_AND_GROWTH_LAW_IN_CELLML_SUBTYPE=44
+  INTEGER(INTG), PARAMETER :: EQUATIONS_SET_MR_AND_EXPONENTIAL_SUBTYPE=45
+  INTEGER(INTG), PARAMETER :: EQUATIONS_SET_REFERENCE_STATE_MR_AND_EXPONENTIAL_SUBTYPE = 46
   
   !Fluid mechanics class
   !  Stokes equations
diff --git a/src/finite_elasticity_routines.F90 b/src/finite_elasticity_routines.F90
index c5b8c450..b16a7afe 100644
--- a/src/finite_elasticity_routines.F90
+++ b/src/finite_elasticity_routines.F90
@@ -669,6 +669,7 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_ELASTICITY_TENSOR(EQUATIONS_SET,DEPENDENT_INT
     REAL(DP) :: TEMPTERM1,TEMPTERM2,VALUE
     REAL(DP), POINTER :: C(:) !Parameters for constitutive laws
     REAL(DP) :: B(6),E(6),DQ_DE(6),Q
+    REAL(DP) :: ONETHIRD_TRACE
     REAL(DP) :: I3EE(6,6) !<Derivative of I3 wrt E
     REAL(DP) :: ADJCC(6,6) !<Derivative of adj(C) wrt C
     REAL(DP) :: AZUE(6,6) !<Derivative of C^-1 wrt E
@@ -751,6 +752,8 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_ELASTICITY_TENSOR(EQUATIONS_SET,DEPENDENT_INT
       STRESS_TENSOR(4)=TEMPTERM1*AZL(2,1)
       STRESS_TENSOR(5)=TEMPTERM1*AZL(3,1)
       STRESS_TENSOR(6)=TEMPTERM1*AZL(3,2)
+
+
       IF(EQUATIONS_SET%specification(3)==EQUATIONS_SET_MOONEY_RIVLIN_ACTIVECONTRACTION_SUBTYPE) THEN
         
         !add active contraction stress values
@@ -779,8 +782,7 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_ELASTICITY_TENSOR(EQUATIONS_SET,DEPENDENT_INT
       ELASTICITY_TENSOR(4,4)=TEMPTERM2
       ELASTICITY_TENSOR(5,5)=TEMPTERM2
       ELASTICITY_TENSOR(6,6)=TEMPTERM2
-      !Add volumetric part of elasticity tensor - p*d(C^-1)/dE.
-      ELASTICITY_TENSOR=ELASTICITY_TENSOR + P*AZUE
+
 
       !Hydrostatic portion of the elasticity tensor (dS/dp)
       HYDRO_ELASTICITY_VOIGT = AZUv
@@ -790,9 +792,56 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_ELASTICITY_TENSOR(EQUATIONS_SET,DEPENDENT_INT
       CALL FINITE_ELASTICITY_PUSH_STRESS_TENSOR(HYDRO_ELASTICITY_VOIGT,DZDNU,Jznu,ERR,ERROR,*999)
       CALL FINITE_ELASTICITY_PUSH_ELASTICITY_TENSOR(ELASTICITY_TENSOR,DZDNU,Jznu,ERR,ERROR,*999)
 
+      !Add volumetric part of elasticity tensor - p*d(C^-1)/dE.
+      ELASTICITY_TENSOR=ELASTICITY_TENSOR + P*AZUE
+
+      ONETHIRD_TRACE=SUM(STRESS_TENSOR(1:3))/3.0_DP
+      STRESS_TENSOR(1:3)=STRESS_TENSOR(1:3)-ONETHIRD_TRACE+P
+
+    CASE(EQUATIONS_SET_ISOTROPIC_EXPONENTIAL_SUBTYPE, EQUATIONS_SET_REFERENCE_STATE_ISO_EXPONENTIAL_SUBTYPE)
+
+      PRESSURE_COMPONENT=DEPENDENT_INTERPOLATED_POINT%INTERPOLATION_PARAMETERS%FIELD_VARIABLE%NUMBER_OF_COMPONENTS
+      P=DEPENDENT_INTERPOLATED_POINT%VALUES(PRESSURE_COMPONENT,NO_PART_DERIV)
+
+      !Form of constitutive model is:
+      ! W=c1/2 (e^(c2*(I1-3)) - 1)
+      ! S = 2*dW/dC
+      I1 = AZL(1,1) + AZL(2,2) + AZL(3,3)
+      TEMPTERM1 = (C(1)*C(2)*EXP(C(2)*(I1-3.0_DP)))*2.0_DP
+      STRESS_TENSOR(1)=TEMPTERM1
+      STRESS_TENSOR(2)=TEMPTERM1
+      STRESS_TENSOR(3)=TEMPTERM1
+      STRESS_TENSOR(4)=0
+      STRESS_TENSOR(5)=0
+      STRESS_TENSOR(6)=0
+
       ! Add volumetric parts.
       STRESS_TENSOR(1:3)=STRESS_TENSOR(1:3)+P
 
+      TEMPTERM1=2.0_DP*TEMPTERM1*C(2)
+      ELASTICITY_TENSOR(1,1)= TEMPTERM1
+      ELASTICITY_TENSOR(2,2)= TEMPTERM1
+      ELASTICITY_TENSOR(3,3)= TEMPTERM1
+      ELASTICITY_TENSOR(4,4)= TEMPTERM1
+      ELASTICITY_TENSOR(5,5)= TEMPTERM1
+      ELASTICITY_TENSOR(6,6)= TEMPTERM1
+
+
+      !Hydrostatic portion of the elasticity tensor (dS/dp)
+      HYDRO_ELASTICITY_VOIGT = AZUv
+
+      ! Do push-forward of 2nd Piola tensor and the material elasticity tensor.
+      CALL FINITE_ELASTICITY_PUSH_STRESS_TENSOR(STRESS_TENSOR,DZDNU,Jznu,ERR,ERROR,*999)
+      CALL FINITE_ELASTICITY_PUSH_STRESS_TENSOR(HYDRO_ELASTICITY_VOIGT,DZDNU,Jznu,ERR,ERROR,*999)
+      CALL FINITE_ELASTICITY_PUSH_ELASTICITY_TENSOR(ELASTICITY_TENSOR,DZDNU,Jznu,ERR,ERROR,*999)
+
+      !Add volumetric part of elasticity tensor - p*d(C^-1)/dE.
+      ELASTICITY_TENSOR=ELASTICITY_TENSOR + P*AZUE
+
+      ONETHIRD_TRACE=SUM(STRESS_TENSOR(1:3))/3.0_DP
+      STRESS_TENSOR(1:3)=STRESS_TENSOR(1:3)-ONETHIRD_TRACE+P
+
+
     CASE(EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE,EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, &
       & EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE)
       PRESSURE_COMPONENT=DEPENDENT_INTERPOLATED_POINT%INTERPOLATION_PARAMETERS%FIELD_VARIABLE%NUMBER_OF_COMPONENTS
@@ -929,7 +978,10 @@ SUBROUTINE FiniteElasticity_FiniteElementJacobianEvaluate(EQUATIONS_SET,ELEMENT_
         nonlinearMatrices=>vectorMatrices%nonlinearMatrices
         jacobianMatrix=>nonlinearMatrices%jacobians(1)%ptr
         IF(jacobianMatrix%updateJacobian) THEN
-          IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) THEN
+          IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE .OR. &
+            & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE .OR. &
+            & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_ISO_EXPONENTIAL_SUBTYPE .OR. &
+            & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MR_AND_EXPONENTIAL_SUBTYPE) THEN
             DEPENDENT_FIELD=>equations%interpolation%geometricField
             GEOMETRIC_FIELD=>equations%interpolation%dependentField
           ELSE
@@ -980,7 +1032,10 @@ SUBROUTINE FiniteElasticity_FiniteElementJacobianEvaluate(EQUATIONS_SET,ELEMENT_
           !Point interpolation pointer
           geometricInterpPoint=>equations%interpolation%geometricInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr
           geometricInterpPointMetrics=>equations%interpolation%geometricInterpPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr
-          IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) THEN
+          IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE .OR. &
+            & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE .OR. &
+            & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_ISO_EXPONENTIAL_SUBTYPE.OR. &
+            & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MR_AND_EXPONENTIAL_SUBTYPE) THEN
             dependentInterpPoint=>equations%interpolation%geometricInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr
             dependentInterpPointMetrics=>equations%interpolation%geometricInterpPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr
             geometricInterpPoint=>equations%interpolation%dependentInterpPoint(FIELD_VAR_TYPE)%ptr
@@ -1513,7 +1568,10 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN
         rhsVector=>vectorMatrices%rhsVector
         vectorMapping =>vectorEquations%vectorMapping
 
-        IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) THEN
+        IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE .OR. &
+          & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE .OR. &
+           & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_ISO_EXPONENTIAL_SUBTYPE .OR. &
+           & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MR_AND_EXPONENTIAL_SUBTYPE) THEN
           DEPENDENT_FIELD  =>equations%interpolation%geometricField
           GEOMETRIC_FIELD  =>equations%interpolation%dependentField
         ELSE
@@ -1570,7 +1628,10 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN
         !Grab interpolation parameters
         FIELD_VARIABLE=>EQUATIONS_SET%equations%vectorEquations%vectorMapping%nonlinearMapping%residualVariables(1)%ptr
         FIELD_VAR_TYPE=FIELD_VARIABLE%VARIABLE_TYPE
-        IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) THEN
+        IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE .OR. &
+          & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE .OR. &
+          & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_ISO_EXPONENTIAL_SUBTYPE .OR. &
+          & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MR_AND_EXPONENTIAL_SUBTYPE) THEN
           GEOMETRIC_INTERPOLATION_PARAMETERS=>equations%interpolation%dependentInterpParameters(FIELD_VAR_TYPE)%ptr
           DEPENDENT_INTERPOLATION_PARAMETERS=>equations%interpolation%geometricInterpParameters(FIELD_U_VARIABLE_TYPE)%ptr
         ELSE
@@ -1631,7 +1692,10 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN
 !       END IF
 
         !Point interpolation pointer
-        IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) THEN
+        IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE .OR. &
+          & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE .OR. &
+          & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_ISO_EXPONENTIAL_SUBTYPE .OR. &
+          & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MR_AND_EXPONENTIAL_SUBTYPE) THEN
           DEPENDENT_INTERPOLATED_POINT=>equations%interpolation%geometricInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr
           DEPENDENT_INTERPOLATED_POINT_METRICS=>equations%interpolation%geometricInterpPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr
           GEOMETRIC_INTERPOLATED_POINT=>equations%interpolation%dependentInterpPoint(FIELD_VAR_TYPE)%ptr
@@ -1667,7 +1731,10 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN
         SELECT CASE(EQUATIONS_SET_SUBTYPE)
         ! ---------------------------------------------------------------
         CASE(EQUATIONS_SET_MOONEY_RIVLIN_ACTIVECONTRACTION_SUBTYPE, &
-          & EQUATIONS_SET_MOONEY_RIVLIN_SUBTYPE)
+          & EQUATIONS_SET_MOONEY_RIVLIN_SUBTYPE, EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE, &
+          & EQUATIONS_SET_ISOTROPIC_EXPONENTIAL_SUBTYPE, EQUATIONS_SET_REFERENCE_STATE_ISO_EXPONENTIAL_SUBTYPE, &
+          & EQUATIONS_SET_MR_AND_EXPONENTIAL_SUBTYPE, EQUATIONS_SET_REFERENCE_STATE_MR_AND_EXPONENTIAL_SUBTYPE )
+
           !Loop over gauss points and add residuals
           DO gauss_idx=1,DEPENDENT_NUMBER_OF_GAUSS_POINTS
             !Interpolate dependent, geometric, fibre and materials fields
@@ -1823,7 +1890,6 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN
         ! ---------------------------------------------------------------
         CASE(EQUATIONS_SET_NO_SUBTYPE,EQUATIONS_SET_MEMBRANE_SUBTYPE, &
           & EQUATIONS_SET_STVENANT_KIRCHOFF_ACTIVECONTRACTION_SUBTYPE, &
-          & EQUATIONS_SET_ISOTROPIC_EXPONENTIAL_SUBTYPE, &
           & EQUATIONS_SET_TRANSVERSE_ISOTROPIC_EXPONENTIAL_SUBTYPE, &
           & EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE,EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, &
           & EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE, &
@@ -2138,7 +2204,10 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN
             
             !Hydrostatic pressure component (skip for membrane problems)
             IF (EQUATIONS_SET_SUBTYPE /= EQUATIONS_SET_MEMBRANE_SUBTYPE) THEN
-              IF(EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) THEN
+              IF(EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE .OR. &
+                & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE .OR. &
+                & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_ISO_EXPONENTIAL_SUBTYPE .OR. &
+                & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MR_AND_EXPONENTIAL_SUBTYPE) THEN
                 HYDROSTATIC_PRESSURE_COMPONENT=GEOMETRIC_FIELD%VARIABLES(var1)%NUMBER_OF_COMPONENTS
                 DEPENDENT_COMPONENT_INTERPOLATION_TYPE=GEOMETRIC_FIELD%VARIABLES(var1)%COMPONENTS( &
                   & HYDROSTATIC_PRESSURE_COMPONENT)%INTERPOLATION_TYPE
@@ -2153,7 +2222,10 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN
                 TEMPTERM1=GAUSS_WEIGHT*Jxxi*(Jznu-Jg)
               ENDIF            
               IF(DEPENDENT_COMPONENT_INTERPOLATION_TYPE==FIELD_NODE_BASED_INTERPOLATION) THEN !node based
-                IF(EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) THEN
+                IF(EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE .OR. &
+                  & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE .OR. &
+                  & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_ISO_EXPONENTIAL_SUBTYPE .OR. &
+                  & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MR_AND_EXPONENTIAL_SUBTYPE) THEN
                   COMPONENT_BASIS=>GEOMETRIC_FIELD%VARIABLES(var1)%COMPONENTS(HYDROSTATIC_PRESSURE_COMPONENT)%DOMAIN% &
                     & TOPOLOGY%ELEMENTS%ELEMENTS(elementNumber)%BASIS
                 ELSE
@@ -2524,7 +2596,10 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN
             
             !Hydrostatic pressure component (skip for membrane problems)
             IF (EQUATIONS_SET_SUBTYPE /= EQUATIONS_SET_MEMBRANE_SUBTYPE) THEN
-              IF(EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) THEN
+              IF(EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE .OR. &
+                & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE .OR. &
+                & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_ISO_EXPONENTIAL_SUBTYPE .OR. &
+                & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MR_AND_EXPONENTIAL_SUBTYPE) THEN
                 HYDROSTATIC_PRESSURE_COMPONENT=GEOMETRIC_FIELD%VARIABLES(var1)%NUMBER_OF_COMPONENTS
                 DEPENDENT_COMPONENT_INTERPOLATION_TYPE=GEOMETRIC_FIELD%VARIABLES(var1)%COMPONENTS( &
                   & HYDROSTATIC_PRESSURE_COMPONENT)%INTERPOLATION_TYPE
@@ -2539,7 +2614,10 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN
                 TEMPTERM1=GAUSS_WEIGHT*Jxxi*(Je-1.0_DP)
               ENDIF            
               IF(DEPENDENT_COMPONENT_INTERPOLATION_TYPE==FIELD_NODE_BASED_INTERPOLATION) THEN !node based
-                IF(EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) THEN
+                IF(EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE .OR. &
+                  & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE .OR. &
+                  & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_ISO_EXPONENTIAL_SUBTYPE .OR. &
+                  & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MR_AND_EXPONENTIAL_SUBTYPE) THEN
                   COMPONENT_BASIS=>GEOMETRIC_FIELD%VARIABLES(var1)%COMPONENTS(HYDROSTATIC_PRESSURE_COMPONENT)%DOMAIN% &
                     & TOPOLOGY%ELEMENTS%ELEMENTS(elementNumber)%BASIS
                 ELSE
@@ -3539,10 +3617,12 @@ SUBROUTINE FiniteElasticity_FiniteElementPreResidualEvaluate(equationsSet,err,er
       CALL Field_VariableGet(dependentField,FIELD_U1_VARIABLE_TYPE,dependentVariable,err,error,*999)
       CALL FiniteElasticity_StressStrainCalculate(equationsSet,EQUATIONS_SET_R_CAUCHY_GREEN_DEFORMATION_TENSOR, &
         & dependentVariable,err,error,*999)
+
     CASE(EQUATIONS_SET_MEMBRANE_SUBTYPE,EQUATIONS_SET_MOONEY_RIVLIN_SUBTYPE, &
       & EQUATIONS_SET_MOONEY_RIVLIN_ACTIVECONTRACTION_SUBTYPE, &
-      & EQUATIONS_SET_STVENANT_KIRCHOFF_ACTIVECONTRACTION_SUBTYPE, &
+      & EQUATIONS_SET_STVENANT_KIRCHOFF_ACTIVECONTRACTION_SUBTYPE,&
       & EQUATIONS_SET_ISOTROPIC_EXPONENTIAL_SUBTYPE, &
+      & EQUATIONS_SET_MR_AND_EXPONENTIAL_SUBTYPE, EQUATIONS_SET_REFERENCE_STATE_MR_AND_EXPONENTIAL_SUBTYPE,&
       & EQUATIONS_SET_TRANSVERSE_ISOTROPIC_EXPONENTIAL_SUBTYPE,EQUATIONS_SET_TRANSVERSE_ISOTROPIC_POLYNOMIAL_SUBTYPE, &
       & EQUATIONS_SET_ANISOTROPIC_POLYNOMIAL_SUBTYPE,EQUATIONS_SET_ANISOTROPIC_POLYNOMIAL_ACTIVE_SUBTYPE, &
       & EQUATIONS_SET_TRANS_ISOTROPIC_ACTIVE_TRANSITION_SUBTYPE,EQUATIONS_SET_TRANSVERSE_ISOTROPIC_ACTIVE_SUBTYPE, &
@@ -3555,6 +3635,8 @@ SUBROUTINE FiniteElasticity_FiniteElementPreResidualEvaluate(equationsSet,err,er
       & EQUATIONS_SET_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE,EQUATIONS_SET_NEARLY_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE, &
       & EQUATIONS_SET_INCOMPRESSIBLE_ELAST_MULTI_COMP_DARCY_SUBTYPE,EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE, &
       & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE, &
+      & EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE, &
+      & EQUATIONS_SET_REFERENCE_STATE_ISO_EXPONENTIAL_SUBTYPE, &
       & EQUATIONS_SET_ELASTICITY_FLUID_PRESSURE_STATIC_INRIA_SUBTYPE, &
       & EQUATIONS_SET_ELASTICITY_FLUID_PRESSURE_HOLMES_MOW_SUBTYPE, &
       & EQUATIONS_SET_ELASTICITY_FLUID_PRES_HOLMES_MOW_ACTIVE_SUBTYPE, &
@@ -3610,6 +3692,8 @@ SUBROUTINE FiniteElasticity_FiniteElementPostResidualEvaluate(EQUATIONS_SET,err,
         & EQUATIONS_SET_MOONEY_RIVLIN_ACTIVECONTRACTION_SUBTYPE, & 
         & EQUATIONS_SET_STVENANT_KIRCHOFF_ACTIVECONTRACTION_SUBTYPE, &
         & EQUATIONS_SET_ISOTROPIC_EXPONENTIAL_SUBTYPE, &
+        & EQUATIONS_SET_MR_AND_EXPONENTIAL_SUBTYPE, &
+        & EQUATIONS_SET_REFERENCE_STATE_MR_AND_EXPONENTIAL_SUBTYPE, &
         & EQUATIONS_SET_TRANSVERSE_ISOTROPIC_EXPONENTIAL_SUBTYPE,EQUATIONS_SET_STANDARD_MONODOMAIN_ELASTICITY_SUBTYPE, &
         & EQUATIONS_SET_ORTHOTROPIC_MATERIAL_COSTA_SUBTYPE,EQUATIONS_SET_COMPRESSIBLE_FINITE_ELASTICITY_SUBTYPE,&
         & EQUATIONS_SET_COMPRESSIBLE_ACTIVECONTRACTION_SUBTYPE,EQUATIONS_SET_TRANSVERSE_ISOTROPIC_ACTIVE_SUBTYPE, &
@@ -3623,6 +3707,8 @@ SUBROUTINE FiniteElasticity_FiniteElementPostResidualEvaluate(EQUATIONS_SET,err,
         & EQUATIONS_SET_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE,EQUATIONS_SET_NEARLY_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE, & 
         & EQUATIONS_SET_INCOMPRESSIBLE_ELAST_MULTI_COMP_DARCY_SUBTYPE,EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE, &
         & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE,  EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE,&
+        & EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE, &
+        & EQUATIONS_SET_REFERENCE_STATE_ISO_EXPONENTIAL_SUBTYPE, &
         & EQUATIONS_SET_ELASTICITY_FLUID_PRESSURE_STATIC_INRIA_SUBTYPE, &
         & EQUATIONS_SET_TRANSVERSE_ISOTROPIC_HUMPHREY_YIN_SUBTYPE,&
         & EQUATIONS_SET_ELASTICITY_FLUID_PRESSURE_HOLMES_MOW_SUBTYPE, &
@@ -4683,8 +4769,13 @@ SUBROUTINE FiniteElasticity_TensorInterpolateXi(equationsSet,tensorEvaluateType,
     !Calculate F=dZ/dNU, the deformation gradient tensor at the xi location
     numberOfDimensions=equationsSet%region%coordinate_system%number_of_dimensions
     numberOfXi=elementBasis%number_of_xi
-    CALL FiniteElasticity_GaussDeformationGradientTensor(dependentInterpolatedPointMetrics, &
+    IF(ASSOCIATED(equations%interpolation%fibreInterpPoint)) THEN
+      CALL FiniteElasticity_GaussDeformationGradientTensor(dependentInterpolatedPointMetrics, &
       & geometricInterpolatedPointMetrics,fibreInterpolatedPoint,dZdNu,err,error,*999)
+    ELSE
+      CALL FiniteElasticity_GaussDeformationGradientTensor(dependentInterpolatedPointMetrics, &
+        & geometricInterpolatedPointMetrics,geometricInterpolatedPoint,dZdNu,err,error,*999)
+    END IF
 
     IF(tensorEvaluateType==EQUATIONS_SET_R_CAUCHY_GREEN_DEFORMATION_TENSOR .OR. &
       & tensorEvaluateType==EQUATIONS_SET_GREEN_LAGRANGE_STRAIN_TENSOR) THEN
@@ -5032,7 +5123,10 @@ SUBROUTINE FiniteElasticity_SurfacePressureResidualEvaluate(EQUATIONS_SET,ELEMEN
     nonlinearMatrices=>vectorEquations%vectorMatrices%nonlinearMatrices
     EQUATIONS_SET_SUBTYPE = EQUATIONS_SET%SPECIFICATION(3)
 
-    IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) THEN
+    IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE .OR. &
+      & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE .OR. &
+      & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_ISO_EXPONENTIAL_SUBTYPE .OR. &
+      & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MR_AND_EXPONENTIAL_SUBTYPE) THEN
       DEPENDENT_FIELD=>equations%interpolation%geometricField
     ELSE
       DEPENDENT_FIELD=>equations%interpolation%dependentField
@@ -5071,7 +5165,10 @@ SUBROUTINE FiniteElasticity_SurfacePressureResidualEvaluate(EQUATIONS_SET,ELEMEN
           DEPENDENT_FACE_BASIS=>DECOMPOSITION%DOMAIN(MESH_COMPONENT_NUMBER)%ptr%TOPOLOGY%FACES%FACES(face_number)%BASIS
           FACE_QUADRATURE_SCHEME=>DEPENDENT_FACE_BASIS%QUADRATURE%QUADRATURE_SCHEME_MAP(BASIS_DEFAULT_QUADRATURE_SCHEME)%ptr
           
-          IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) THEN
+          IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE .OR. &
+            & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE .OR. &
+            & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_ISO_EXPONENTIAL_SUBTYPE .OR. &
+            & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MR_AND_EXPONENTIAL_SUBTYPE) THEN
             FACE_DEPENDENT_INTERPOLATION_PARAMETERS=>equations%interpolation%geometricInterpParameters(FIELD_VAR_U_TYPE)%ptr
             FACE_DEPENDENT_INTERPOLATED_POINT=>equations%interpolation%geometricInterpPoint(FIELD_VAR_U_TYPE)%ptr
             FACE_DEPENDENT_INTERPOLATED_POINT_METRICS=>equations%interpolation% &
@@ -6038,7 +6135,19 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_CAUCHY_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPO
       & EQUATIONS_SET_INCOMPRESSIBLE_FINITE_ELASTICITY_DARCY_SUBTYPE,EQUATIONS_SET_STANDARD_MONODOMAIN_ELASTICITY_SUBTYPE, &
       & EQUATIONS_SET_1D3D_MONODOMAIN_ELASTICITY_SUBTYPE,EQUATIONS_SET_TRANSVERSE_ISOTROPIC_POLYNOMIAL_SUBTYPE, &
       & EQUATIONS_SET_MONODOMAIN_ELASTICITY_W_TITIN_SUBTYPE,EQUATIONS_SET_MONODOMAIN_ELASTICITY_VELOCITY_SUBTYPE, &
-      & EQUATIONS_SET_TRANSVERSE_ISOTROPIC_ACTIVE_SUBTYPE,EQUATIONS_SET_MR_AND_GROWTH_LAW_IN_CELLML_SUBTYPE)
+      & EQUATIONS_SET_TRANSVERSE_ISOTROPIC_ACTIVE_SUBTYPE,EQUATIONS_SET_MR_AND_GROWTH_LAW_IN_CELLML_SUBTYPE, &
+      & EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE, &
+      & EQUATIONS_SET_MR_AND_EXPONENTIAL_SUBTYPE, EQUATIONS_SET_REFERENCE_STATE_MR_AND_EXPONENTIAL_SUBTYPE)
+
+
+      IF(EQUATIONS_SET%specification(3)==EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE) THEN
+        PRESSURE_COMPONENT=GEOMETRIC_INTERPOLATED_POINT%INTERPOLATION_PARAMETERS%FIELD_VARIABLE%NUMBER_OF_COMPONENTS
+        P=GEOMETRIC_INTERPOLATED_POINT%VALUES(PRESSURE_COMPONENT,NO_PART_DERIV)
+      ELSE
+        PRESSURE_COMPONENT=DEPENDENT_INTERPOLATED_POINT%INTERPOLATION_PARAMETERS%FIELD_VARIABLE%NUMBER_OF_COMPONENTS
+        P=DEPENDENT_INTERPOLATED_POINT%VALUES(PRESSURE_COMPONENT,NO_PART_DERIV)
+      ENDIF
+
       !Form of constitutive model is:
       ! W=c1*(I1-3)+c2*(I2-3)+p*(I3-1)
       !Also assumed I3 = det(AZL) = 1.0
@@ -6367,11 +6476,6 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_CAUCHY_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPO
         PIOLA_TENSOR(1,1)=PIOLA_TENSOR(1,1)+VALUE
       ENDIF
 
-    CASE(EQUATIONS_SET_ISOTROPIC_EXPONENTIAL_SUBTYPE)
-      !Form of constitutive model is:
-      ! W=c1/2 (e^(c2*(I1-3)) - 1)
-      ! S = 2*dW/dC + 2pC^-1
-      PIOLA_TENSOR=C(1)*C(2)*EXP(C(2)*(AZL(1,1)+AZL(2,2)+AZL(3,3)-3.0_DP))*IDENTITY+2.0_DP*P*AZU
     CASE(EQUATIONS_SET_ELASTICITY_FLUID_PRESSURE_STATIC_INRIA_SUBTYPE)
       !C(1)=Mooney Rivlin parameter
       !C(2)=Mooney Rivlin parameter
@@ -7032,7 +7136,7 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPO
       !Calculate isochoric fictitious 2nd Piola tensor (in Voigt form)
       I1=AZL(1,1)+AZL(2,2)+AZL(3,3)
       TEMPTERM1=-2.0_DP*C(2)
-      TEMPTERM2=2.0_DP*(C(1)+I1*C(2))
+      TEMPTERM2= 2.0_DP*(C(1)+I1*C(2))
       STRESS_TENSOR(1)=TEMPTERM1*AZL(1,1)+TEMPTERM2
       STRESS_TENSOR(2)=TEMPTERM1*AZL(2,2)+TEMPTERM2
       STRESS_TENSOR(3)=TEMPTERM1*AZL(3,3)+TEMPTERM2
@@ -7060,6 +7164,75 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPO
       ONETHIRD_TRACE=SUM(STRESS_TENSOR(1:3))/3.0_DP
       STRESS_TENSOR(1:3)=STRESS_TENSOR(1:3)-ONETHIRD_TRACE+P
 
+
+    CASE(EQUATIONS_SET_ISOTROPIC_EXPONENTIAL_SUBTYPE, EQUATIONS_SET_REFERENCE_STATE_ISO_EXPONENTIAL_SUBTYPE)
+
+      IF(EQUATIONS_SET%specification(3)==EQUATIONS_SET_REFERENCE_STATE_ISO_EXPONENTIAL_SUBTYPE) THEN
+        PRESSURE_COMPONENT=GEOMETRIC_INTERPOLATED_POINT%INTERPOLATION_PARAMETERS%FIELD_VARIABLE%NUMBER_OF_COMPONENTS
+        P=GEOMETRIC_INTERPOLATED_POINT%VALUES(PRESSURE_COMPONENT,NO_PART_DERIV)
+      ELSE
+        PRESSURE_COMPONENT=DEPENDENT_INTERPOLATED_POINT%INTERPOLATION_PARAMETERS%FIELD_VARIABLE%NUMBER_OF_COMPONENTS
+        P=DEPENDENT_INTERPOLATED_POINT%VALUES(PRESSURE_COMPONENT,NO_PART_DERIV)
+      ENDIF
+      !Form of constitutive model is:
+      ! W=c1/2 (e^(c2*(I1-3)) - 1)
+      ! S = 2*dW/dC + 2pC^-1
+      I1 = AZL(1,1) + AZL(2,2) + AZL(3,3)
+      TEMPTERM1 = (C(1)*C(2)*EXP(C(2)*(I1-3.0_DP)))*2.0_DP
+      STRESS_TENSOR(1)=TEMPTERM1
+      STRESS_TENSOR(2)=TEMPTERM1
+      STRESS_TENSOR(3)=TEMPTERM1
+      STRESS_TENSOR(4)=0
+      STRESS_TENSOR(5)=0
+      STRESS_TENSOR(6)=0
+
+
+      CALL FINITE_ELASTICITY_PUSH_STRESS_TENSOR(STRESS_TENSOR,MOD_DZDNU,Jznu,err,error,*999)
+      !Calculate isochoric Cauchy tensor (the deviatoric part) and add the volumetric part (the hydrostatic pressure).
+      ONETHIRD_TRACE=SUM(STRESS_TENSOR(1:3))/3.0_DP
+      STRESS_TENSOR(1:3)=STRESS_TENSOR(1:3)-ONETHIRD_TRACE+P
+
+    CASE(EQUATIONS_SET_MR_AND_EXPONENTIAL_SUBTYPE, EQUATIONS_SET_REFERENCE_STATE_MR_AND_EXPONENTIAL_SUBTYPE)
+
+      IF(EQUATIONS_SET%specification(3)== EQUATIONS_SET_REFERENCE_STATE_MR_AND_EXPONENTIAL_SUBTYPE) THEN
+        PRESSURE_COMPONENT=GEOMETRIC_INTERPOLATED_POINT%INTERPOLATION_PARAMETERS%FIELD_VARIABLE%NUMBER_OF_COMPONENTS
+        P=GEOMETRIC_INTERPOLATED_POINT%VALUES(PRESSURE_COMPONENT,NO_PART_DERIV)
+      ELSE
+        PRESSURE_COMPONENT=DEPENDENT_INTERPOLATED_POINT%INTERPOLATION_PARAMETERS%FIELD_VARIABLE%NUMBER_OF_COMPONENTS
+        P=DEPENDENT_INTERPOLATED_POINT%VALUES(PRESSURE_COMPONENT,NO_PART_DERIV)
+      ENDIF
+
+
+      IF (NINT(C(5)) .EQ. 1) THEN
+        I1=AZL(1,1)+AZL(2,2)+AZL(3,3)
+        TEMPTERM1=-2.0_DP*C(2)
+        TEMPTERM2= 2.0_DP*(C(1)+I1*C(2))
+        STRESS_TENSOR(1)=TEMPTERM1*AZL(1,1)+TEMPTERM2
+        STRESS_TENSOR(2)=TEMPTERM1*AZL(2,2)+TEMPTERM2
+        STRESS_TENSOR(3)=TEMPTERM1*AZL(3,3)+TEMPTERM2
+        STRESS_TENSOR(4)=TEMPTERM1*AZL(2,1)
+        STRESS_TENSOR(5)=TEMPTERM1*AZL(3,1)
+        STRESS_TENSOR(6)=TEMPTERM1*AZL(3,2)
+
+      ELSE
+        !Form of constitutive model is:
+        ! W=c1/2 (e^(c2*(I1-3)) - 1)
+        ! S = 2*dW/dC + 2pC^-1
+        I1 = AZL(1,1) + AZL(2,2) + AZL(3,3)
+        TEMPTERM1 = (C(3)*C(4)*EXP(C(4)*(I1-3.0_DP)))*2.0_DP
+        STRESS_TENSOR(1)=TEMPTERM1
+        STRESS_TENSOR(2)=TEMPTERM1
+        STRESS_TENSOR(3)=TEMPTERM1
+        STRESS_TENSOR(4)=0
+        STRESS_TENSOR(5)=0
+        STRESS_TENSOR(6)=0
+      END IF
+
+      CALL FINITE_ELASTICITY_PUSH_STRESS_TENSOR(STRESS_TENSOR,MOD_DZDNU,Jznu,err,error,*999)
+      !Calculate isochoric Cauchy tensor (the deviatoric part) and add the volumetric part (the hydrostatic pressure).
+      ONETHIRD_TRACE=SUM(STRESS_TENSOR(1:3))/3.0_DP
+      STRESS_TENSOR(1:3)=STRESS_TENSOR(1:3)-ONETHIRD_TRACE+P
+
     CASE(EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE,EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, &
       & EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE)
       IF(EQUATIONS_SET%specification(3)==EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) THEN
@@ -7386,6 +7559,7 @@ SUBROUTINE FINITE_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET
         & EQUATIONS_SET_STVENANT_KIRCHOFF_ACTIVECONTRACTION_SUBTYPE, &
         & EQUATIONS_SET_ISOTROPIC_EXPONENTIAL_SUBTYPE,EQUATIONS_SET_TRANSVERSE_ISOTROPIC_ACTIVE_SUBTYPE, &
         & EQUATIONS_SET_TRANS_ISOTROPIC_ACTIVE_TRANSITION_SUBTYPE, &
+        & EQUATIONS_SET_MR_AND_EXPONENTIAL_SUBTYPE, EQUATIONS_SET_REFERENCE_STATE_MR_AND_EXPONENTIAL_SUBTYPE, &
         & EQUATIONS_SET_TRANSVERSE_ISOTROPIC_EXPONENTIAL_SUBTYPE,EQUATIONS_SET_TRANSVERSE_ISOTROPIC_POLYNOMIAL_SUBTYPE, &
         & EQUATIONS_SET_ANISOTROPIC_POLYNOMIAL_SUBTYPE,EQUATIONS_SET_ANISOTROPIC_POLYNOMIAL_ACTIVE_SUBTYPE, &
         & EQUATIONS_SET_ORTHOTROPIC_MATERIAL_COSTA_SUBTYPE,EQUATIONS_SET_COMPRESSIBLE_FINITE_ELASTICITY_SUBTYPE, &
@@ -7397,6 +7571,8 @@ SUBROUTINE FINITE_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET
         & EQUATIONS_SET_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE,EQUATIONS_SET_NEARLY_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE, & 
         & EQUATIONS_SET_INCOMPRESSIBLE_ELAST_MULTI_COMP_DARCY_SUBTYPE,EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE, &
         & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE, &
+        & EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE, &
+        & EQUATIONS_SET_REFERENCE_STATE_ISO_EXPONENTIAL_SUBTYPE, &
         & EQUATIONS_SET_CONSTITUTIVE_LAW_IN_CELLML_EVALUATE_SUBTYPE, &
         & EQUATIONS_SET_CONSTITUTIVE_AND_GROWTH_LAW_IN_CELLML_SUBTYPE, &
         & EQUATIONS_SET_GROWTH_LAW_IN_CELLML_SUBTYPE, &
@@ -7409,6 +7585,7 @@ SUBROUTINE FINITE_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET
         & EQUATIONS_SET_HOLZAPFEL_OGDEN_ACTIVECONTRACTION_SUBTYPE, &
         & EQUATIONS_SET_RATE_BASED_SMOOTH_MODEL_SUBTYPE,EQUATIONS_SET_COMPRESSIBLE_RATE_BASED_SMOOTH_MODEL_SUBTYPE, &
         & EQUATIONS_SET_RATE_BASED_GROWTH_MODEL_SUBTYPE,EQUATIONS_SET_COMPRESSIBLE_RATE_BASED_GROWTH_MODEL_SUBTYPE)
+
         SELECT CASE(EQUATIONS_SET_SETUP%SETUP_TYPE)
         CASE(EQUATIONS_SET_SETUP_INITIAL_TYPE)
           SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE)
@@ -7481,8 +7658,12 @@ SUBROUTINE FINITE_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET
             SELECT CASE(EQUATIONS_SET_SUBTYPE)
             CASE(EQUATIONS_SET_MOONEY_RIVLIN_SUBTYPE, &
               & EQUATIONS_SET_MOONEY_RIVLIN_ACTIVECONTRACTION_SUBTYPE, &
+              & EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE, &
+              & EQUATIONS_SET_REFERENCE_STATE_ISO_EXPONENTIAL_SUBTYPE, &
               & EQUATIONS_SET_STVENANT_KIRCHOFF_ACTIVECONTRACTION_SUBTYPE, &
               & EQUATIONS_SET_ISOTROPIC_EXPONENTIAL_SUBTYPE, &
+              & EQUATIONS_SET_MR_AND_EXPONENTIAL_SUBTYPE, &
+              & EQUATIONS_SET_REFERENCE_STATE_MR_AND_EXPONENTIAL_SUBTYPE, &
               & EQUATIONS_SET_COMPRESSIBLE_FINITE_ELASTICITY_SUBTYPE,&
               & EQUATIONS_SET_COMPRESSIBLE_ACTIVECONTRACTION_SUBTYPE,&
               & EQUATIONS_SET_NO_SUBTYPE, &
@@ -7492,6 +7673,7 @@ SUBROUTINE FINITE_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET
               & EQUATIONS_SET_ELASTICITY_FLUID_PRES_HOLMES_MOW_ACTIVE_SUBTYPE, &
               & EQUATIONS_SET_NEARLY_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE, & 
               & EQUATIONS_SET_MR_AND_GROWTH_LAW_IN_CELLML_SUBTYPE)
+
               ! pass, fibre field isn't required as the constitutive relation is isotropic
             CASE(EQUATIONS_SET_ORTHOTROPIC_MATERIAL_COSTA_SUBTYPE, &
               & EQUATIONS_SET_TRANSVERSE_ISOTROPIC_EXPONENTIAL_SUBTYPE, &
@@ -7591,12 +7773,17 @@ SUBROUTINE FINITE_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET
             & EQUATIONS_SET_ACTIVECONTRACTION_SUBTYPE, EQUATIONS_SET_NO_SUBTYPE,EQUATIONS_SET_MEMBRANE_SUBTYPE, &
             & EQUATIONS_SET_ORTHOTROPIC_MATERIAL_HOLZAPFEL_OGDEN_SUBTYPE,EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE, &
             & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE, &
+            & EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE, &
+            & EQUATIONS_SET_REFERENCE_STATE_ISO_EXPONENTIAL_SUBTYPE, &
             & EQUATIONS_SET_TRANSVERSE_ISOTROPIC_HUMPHREY_YIN_SUBTYPE, &
             & EQUATIONS_SET_STANDARD_MONODOMAIN_ELASTICITY_SUBTYPE,EQUATIONS_SET_1D3D_MONODOMAIN_ELASTICITY_SUBTYPE, &
             & EQUATIONS_SET_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE,EQUATIONS_SET_NEARLY_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE, &
             & EQUATIONS_SET_MONODOMAIN_ELASTICITY_W_TITIN_SUBTYPE, &
             & EQUATIONS_SET_MONODOMAIN_ELASTICITY_VELOCITY_SUBTYPE, &
-            & EQUATIONS_SET_HOLZAPFEL_OGDEN_ACTIVECONTRACTION_SUBTYPE)
+            & EQUATIONS_SET_HOLZAPFEL_OGDEN_ACTIVECONTRACTION_SUBTYPE, &
+            & EQUATIONS_SET_MR_AND_EXPONENTIAL_SUBTYPE, &
+            & EQUATIONS_SET_REFERENCE_STATE_MR_AND_EXPONENTIAL_SUBTYPE)
+
             SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE)
             CASE(EQUATIONS_SET_SETUP_START_ACTION)
               IF(EQUATIONS_SET%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN
@@ -10167,11 +10354,17 @@ SUBROUTINE FINITE_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET
               NUMBER_OF_FLUID_COMPONENTS=0
               SELECT CASE(EQUATIONS_SET_SUBTYPE)
               CASE(EQUATIONS_SET_MOONEY_RIVLIN_SUBTYPE,EQUATIONS_SET_NO_SUBTYPE, &
+                & EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE, &
+                & EQUATIONS_SET_ISOTROPIC_EXPONENTIAL_SUBTYPE, &
+                & EQUATIONS_SET_REFERENCE_STATE_ISO_EXPONENTIAL_SUBTYPE, &
                 & EQUATIONS_SET_MOONEY_RIVLIN_ACTIVECONTRACTION_SUBTYPE, &
                 & EQUATIONS_SET_INCOMPRESSIBLE_FINITE_ELASTICITY_DARCY_SUBTYPE, &
                 & EQUATIONS_SET_STANDARD_MONODOMAIN_ELASTICITY_SUBTYPE, &
                 & EQUATIONS_SET_MR_AND_GROWTH_LAW_IN_CELLML_SUBTYPE)
                 NUMBER_OF_COMPONENTS=2;
+              CASE(EQUATIONS_SET_MR_AND_EXPONENTIAL_SUBTYPE, &
+                & EQUATIONS_SET_REFERENCE_STATE_MR_AND_EXPONENTIAL_SUBTYPE)
+                NUMBER_OF_COMPONENTS=5;
               CASE(EQUATIONS_SET_ACTIVE_STRAIN_SUBTYPE)
                 NUMBER_OF_COMPONENTS=8;
               CASE(EQUATIONS_SET_MULTISCALE_ACTIVE_STRAIN_SUBTYPE)
@@ -10197,8 +10390,6 @@ SUBROUTINE FINITE_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET
                 ENDIF
               CASE(EQUATIONS_SET_STVENANT_KIRCHOFF_ACTIVECONTRACTION_SUBTYPE)
                 NUMBER_OF_COMPONENTS=2;
-              CASE(EQUATIONS_SET_ISOTROPIC_EXPONENTIAL_SUBTYPE)
-                NUMBER_OF_COMPONENTS=2;
               CASE(EQUATIONS_SET_TRANSVERSE_ISOTROPIC_EXPONENTIAL_SUBTYPE)
                 NUMBER_OF_COMPONENTS=5;
               CASE(EQUATIONS_SET_TRANSVERSE_ISOTROPIC_POLYNOMIAL_SUBTYPE)
@@ -10820,6 +11011,8 @@ SUBROUTINE FiniteElasticity_EquationsSetSolutionMethodSet(EQUATIONS_SET,SOLUTION
           & EQUATIONS_SET_INCOMPRESSIBLE_ELASTICITY_DRIVEN_DARCY_SUBTYPE, &
           & EQUATIONS_SET_INCOMPRESSIBLE_ELAST_MULTI_COMP_DARCY_SUBTYPE,EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE, &
           & EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE, &
+          & EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE, &
+          & EQUATIONS_SET_REFERENCE_STATE_ISO_EXPONENTIAL_SUBTYPE, &
           & EQUATIONS_SET_CONSTITUTIVE_LAW_IN_CELLML_EVALUATE_SUBTYPE, &
           & EQUATIONS_SET_CONSTITUTIVE_AND_GROWTH_LAW_IN_CELLML_SUBTYPE, &
           & EQUATIONS_SET_GROWTH_LAW_IN_CELLML_SUBTYPE, &
@@ -10838,7 +11031,9 @@ SUBROUTINE FiniteElasticity_EquationsSetSolutionMethodSet(EQUATIONS_SET,SOLUTION
           & EQUATIONS_SET_1D3D_MONODOMAIN_ACTIVE_STRAIN_SUBTYPE, &
           & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, &
           & EQUATIONS_SET_RATE_BASED_SMOOTH_MODEL_SUBTYPE,EQUATIONS_SET_COMPRESSIBLE_RATE_BASED_SMOOTH_MODEL_SUBTYPE, &
-          & EQUATIONS_SET_RATE_BASED_GROWTH_MODEL_SUBTYPE,EQUATIONS_SET_COMPRESSIBLE_RATE_BASED_GROWTH_MODEL_SUBTYPE)
+          & EQUATIONS_SET_RATE_BASED_GROWTH_MODEL_SUBTYPE,EQUATIONS_SET_COMPRESSIBLE_RATE_BASED_GROWTH_MODEL_SUBTYPE, &
+          & EQUATIONS_SET_MR_AND_EXPONENTIAL_SUBTYPE, &
+          & EQUATIONS_SET_REFERENCE_STATE_MR_AND_EXPONENTIAL_SUBTYPE)
         SELECT CASE(SOLUTION_METHOD)
         CASE(EQUATIONS_SET_FEM_SOLUTION_METHOD)
           EQUATIONS_SET%SOLUTION_METHOD=EQUATIONS_SET_FEM_SOLUTION_METHOD
@@ -10901,7 +11096,8 @@ SUBROUTINE FiniteElasticity_EquationsSetSpecificationSet(equationsSet,specificat
       CASE(EQUATIONS_SET_MEMBRANE_SUBTYPE,EQUATIONS_SET_MOONEY_RIVLIN_SUBTYPE, &
           & EQUATIONS_SET_MOONEY_RIVLIN_ACTIVECONTRACTION_SUBTYPE, &
           & EQUATIONS_SET_STVENANT_KIRCHOFF_ACTIVECONTRACTION_SUBTYPE, &
-          & EQUATIONS_SET_ISOTROPIC_EXPONENTIAL_SUBTYPE,EQUATIONS_SET_TRANSVERSE_ISOTROPIC_ACTIVE_SUBTYPE, &
+          & EQUATIONS_SET_ISOTROPIC_EXPONENTIAL_SUBTYPE,&
+          & EQUATIONS_SET_TRANSVERSE_ISOTROPIC_ACTIVE_SUBTYPE, &
           & EQUATIONS_SET_TRANS_ISOTROPIC_ACTIVE_TRANSITION_SUBTYPE, &
           & EQUATIONS_SET_TRANSVERSE_ISOTROPIC_EXPONENTIAL_SUBTYPE,EQUATIONS_SET_TRANSVERSE_ISOTROPIC_POLYNOMIAL_SUBTYPE, &
           & EQUATIONS_SET_ANISOTROPIC_POLYNOMIAL_SUBTYPE,EQUATIONS_SET_ANISOTROPIC_POLYNOMIAL_ACTIVE_SUBTYPE, &
@@ -10914,6 +11110,8 @@ SUBROUTINE FiniteElasticity_EquationsSetSpecificationSet(equationsSet,specificat
           & EQUATIONS_SET_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE,EQUATIONS_SET_NEARLY_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE, &
           & EQUATIONS_SET_INCOMPRESSIBLE_ELAST_MULTI_COMP_DARCY_SUBTYPE,EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE, &
           & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE, &
+          & EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE, &
+          & EQUATIONS_SET_REFERENCE_STATE_ISO_EXPONENTIAL_SUBTYPE, &
           & EQUATIONS_SET_CONSTITUTIVE_LAW_IN_CELLML_EVALUATE_SUBTYPE, &
           & EQUATIONS_SET_CONSTITUTIVE_AND_GROWTH_LAW_IN_CELLML_SUBTYPE, &
           & EQUATIONS_SET_GROWTH_LAW_IN_CELLML_SUBTYPE, &
@@ -10928,7 +11126,9 @@ SUBROUTINE FiniteElasticity_EquationsSetSpecificationSet(equationsSet,specificat
           & EQUATIONS_SET_1D3D_MONODOMAIN_ACTIVE_STRAIN_SUBTYPE, &
           & EQUATIONS_SET_HOLZAPFEL_OGDEN_ACTIVECONTRACTION_SUBTYPE, &
           & EQUATIONS_SET_RATE_BASED_SMOOTH_MODEL_SUBTYPE,EQUATIONS_SET_COMPRESSIBLE_RATE_BASED_SMOOTH_MODEL_SUBTYPE, &
-          & EQUATIONS_SET_RATE_BASED_GROWTH_MODEL_SUBTYPE,EQUATIONS_SET_COMPRESSIBLE_RATE_BASED_GROWTH_MODEL_SUBTYPE)
+          & EQUATIONS_SET_RATE_BASED_GROWTH_MODEL_SUBTYPE,EQUATIONS_SET_COMPRESSIBLE_RATE_BASED_GROWTH_MODEL_SUBTYPE, &
+          & EQUATIONS_SET_MR_AND_EXPONENTIAL_SUBTYPE, &
+          & EQUATIONS_SET_REFERENCE_STATE_MR_AND_EXPONENTIAL_SUBTYPE)
         !Set full specification
         IF(ALLOCATED(equationsSet%specification)) THEN
           CALL FlagError("Equations set specification is already allocated.",err,error,*999)
diff --git a/src/opencmiss_iron.F90 b/src/opencmiss_iron.F90
index 300cf38f..a5deeceb 100644
--- a/src/opencmiss_iron.F90
+++ b/src/opencmiss_iron.F90
@@ -2518,6 +2518,10 @@ MODULE OpenCMISS_Iron
   INTEGER(INTG), PARAMETER :: CMFE_EQUATIONS_SET_MOONEY_RIVLIN_SUBTYPE = EQUATIONS_SET_MOONEY_RIVLIN_SUBTYPE !< Mooney-Rivlin constitutive law for finite elasticity equations set subtype \see OpenCMISS_EquationsSetSubtypes,OpenCMISS
   INTEGER(INTG), PARAMETER :: CMFE_EQUATIONS_SET_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE = &
   & EQUATIONS_SET_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE !< Incompressible Mooney-Rivlin constitutive law for finite elasticity equations set subtype \see OpenCMISS_EquationsSetSubtypes,OpenCMISS
+  INTEGER(INTG), PARAMETER :: CMFE_EQUATIONS_SET_MR_AND_EXPONENTIAL_SUBTYPE = &
+  & EQUATIONS_SET_MR_AND_EXPONENTIAL_SUBTYPE !< mixing mooney rivlin with exponential constitutive law
+  INTEGER(INTG), PARAMETER :: CMFE_EQUATIONS_SET_REFERENCE_STATE_MR_AND_EXPONENTIAL_SUBTYPE = &
+  & EQUATIONS_SET_REFERENCE_STATE_MR_AND_EXPONENTIAL_SUBTYPE !< mixing mooney rivlin with exponential constitutive law for backward formulation
   INTEGER(INTG), PARAMETER :: CMFE_EQUATIONS_SET_NEARLY_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE = &
   & EQUATIONS_SET_NEARLY_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE !< Nearly Incompressible Mooney-Rivlin constitutive law for finite elasticity equations set subtype \see OpenCMISS_EquationsSetSubtypes,OpenCMISS
   INTEGER(INTG), PARAMETER :: CMFE_EQUATIONS_SET_MOONEY_RIVLIN_ACTIVECONTRACTION_SUBTYPE = &
@@ -2526,7 +2530,8 @@ MODULE OpenCMISS_Iron
     & EQUATIONS_SET_STVENANT_KIRCHOFF_ACTIVECONTRACTION_SUBTYPE !< St Venant Kirchoff constitutive law with steady-state active contraction for finite elasticity equations set subtype \see OpenCMISS_EquationsSetSubtypes,OpenCMISS
   INTEGER(INTG), PARAMETER :: CMFE_EQUATIONS_SET_ACTIVECONTRACTION_SUBTYPE =&
     & EQUATIONS_SET_ACTIVECONTRACTION_SUBTYPE !< Active contraction/costa-based law with quasistatic time loop for finite elasticity equations set subtype \see OpenCMISS_EquationsSetSubtypes,OpenCMISS
-  INTEGER(INTG), PARAMETER :: CMFE_EQUATIONS_SET_ISOTROPIC_EXPONENTIAL_SUBTYPE = EQUATIONS_SET_ISOTROPIC_EXPONENTIAL_SUBTYPE !< Isotropic exponential constitutive law for finite elasticity equations set subtype \see OpenCMISS_EquationsSetSubtypes,OpenCMISS
+  INTEGER(INTG), PARAMETER :: CMFE_EQUATIONS_SET_ISOTROPIC_EXPONENTIAL_SUBTYPE = &
+    & EQUATIONS_SET_ISOTROPIC_EXPONENTIAL_SUBTYPE !< Isotropic exponential constitutive law for finite elasticity equations set subtype \see OpenCMISS_EquationsSetSubtypes,OpenCMISS
   INTEGER(INTG), PARAMETER :: CMFE_EQUATIONS_SET_TRANSVERSE_ISOTROPIC_EXPONENTIAL_SUBTYPE = &
     & EQUATIONS_SET_TRANSVERSE_ISOTROPIC_EXPONENTIAL_SUBTYPE !< Transverse isotropic exponential constitutive law for finite elasticity equations set subtype \see OpenCMISS_EquationsSetSubtypes,OpenCMISS
   INTEGER(INTG), PARAMETER :: CMFE_EQUATIONS_SET_TRANSVERSE_ISOTROPIC_ACTIVE_SUBTYPE = &
@@ -3083,6 +3088,8 @@ MODULE OpenCMISS_Iron
     & CMFE_EQUATIONS_SET_SHELL_SUBTYPE, &
     & CMFE_EQUATIONS_SET_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE,CMFE_EQUATIONS_SET_NEARLY_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE, &
     & CMFE_EQUATIONS_SET_MOONEY_RIVLIN_SUBTYPE, &
+    & CMFE_EQUATIONS_SET_MR_AND_EXPONENTIAL_SUBTYPE,&
+    & CMFE_EQUATIONS_SET_REFERENCE_STATE_MR_AND_EXPONENTIAL_SUBTYPE,&
     & CMFE_EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE, CMFE_EQUATIONS_SET_ISOTROPIC_EXPONENTIAL_SUBTYPE, &
     & CMFE_EQUATIONS_SET_ACTIVECONTRACTION_SUBTYPE,CMFE_EQUATIONS_SET_MOONEY_RIVLIN_ACTIVECONTRACTION_SUBTYPE, &
     & CMFE_EQUATIONS_SET_COMPRESSIBLE_ACTIVECONTRACTION_SUBTYPE,CMFE_EQUATIONS_SET_TRANSVERSE_ISOTROPIC_ACTIVE_SUBTYPE, &