diff --git a/temoa_model/temoa_initialize.py b/temoa_model/temoa_initialize.py
index 0d3d9bd3..3eb6bbea 100644
--- a/temoa_model/temoa_initialize.py
+++ b/temoa_model/temoa_initialize.py
@@ -837,8 +837,9 @@ def RegionalGlobalInitializedIndices ( M ):
 	for n in range(1,len(M.regions)+1):
 		regional_perms = permutations(M.regions,n)
 		for i in regional_perms:
-			indices.add("-".join(i))
+			indices.add("+".join(i))
 	indices.add('global')
+	indices = indices.union(M.RegionalIndices)
 
 	return indices
 
diff --git a/temoa_model/temoa_rules.py b/temoa_model/temoa_rules.py
index e9a9b930..64cd401e 100644
--- a/temoa_model/temoa_rules.py
+++ b/temoa_model/temoa_rules.py
@@ -1531,13 +1531,18 @@ def EmissionLimit_Constraint(M, r, p, e):
 """
     emission_limit = M.EmissionLimit[r, p, e]
 
-     # r can be an individual region (r='US'), or a combination of regions separated by hyphen (r='Mexico-US-Canada'), or 'global'.
-     # Note that regions!=M.regions. We iterate over regions to find actural_emissions and actual_emissions_annual.
-    regions = set(r.split("-"))
+     # r can be an individual region (r='US'), or a combination of regions separated by a + (r='Mexico+US+Canada'), or 'global'.
+     # Note that regions!=M.regions. We iterate over regions to find actual_emissions and actual_emissions_annual.
+    
 
     # if r == 'global', the constraint is system-wide
-    if regions == {'global'}:
+
+    if r == 'global':
       regions = M.regions
+    elif '+' in r:
+      regions = r.split('+')
+    else:
+      regions = [r]
 
 
     actual_emissions = sum(
@@ -1672,17 +1677,19 @@ def MaxActivity_Constraint(M, r, p, t):
    \forall \{r, p, t \in T^{a}\} \in \Theta_{\text{MaxActivity}}
 
 """
-	# r can be an individual region (r='US'), or a combination of regions separated by hyphen (r='Mexico-US-Canada'), or 'global'.
+	# r can be an individual region (r='US'), or a combination of regions separated by a + (r='Mexico+US+Canada'), or 'global'.
     # if r == 'global', the constraint is system-wide
     if r == 'global':
       reg = M.regions
+    elif '+' in r:
+      reg = r.split('+')
     else:
       reg = [r]
 
     try:
       activity_rpt = sum(
           M.V_FlowOut[r, p, s, d, S_i, t, S_v, S_o]
-          for r in reg if '-' not in r
+          for r in reg
           for S_v in M.processVintages[r, p, t]
           for S_i in M.processInputs[r, p, t, S_v]
           for S_o in M.ProcessOutputsByInput[r, p, t, S_v, S_i]
@@ -1692,7 +1699,7 @@ def MaxActivity_Constraint(M, r, p, t):
     except:
       activity_rpt = sum(
           M.V_FlowOutAnnual[r, p, S_i, t, S_v, S_o]
-          for r in reg if '-' not in r
+          for r in reg
           for S_v in M.processVintages[r, p, t]
           for S_i in M.processInputs[r, p, t, S_v]
           for S_o in M.ProcessOutputsByInput[r, p, t, S_v, S_i]
@@ -1725,17 +1732,19 @@ def MinActivity_Constraint(M, r, p, t):
    \forall \{r, p, t \in T^{a}\} \in \Theta_{\text{MinActivity}}
 
 """
-	# r can be an individual region (r='US'), or a combination of regions separated by hyphen (r='Mexico-US-Canada'), or 'global'.
+	# r can be an individual region (r='US'), or a combination of regions separated by a + (r='Mexico+US+Canada'), or 'global'.
     # if r == 'global', the constraint is system-wide
     if r == 'global':
       reg = M.regions
+    elif '+' in r:
+      reg = r.split('+')
     else:
       reg = [r]
 
     try:
       activity_rpt = sum(
           M.V_FlowOut[r, p, s, d, S_i, t, S_v, S_o]
-          for r in reg if '-' not in r
+          for r in reg
           for S_v in M.processVintages[r, p, t]
           for S_i in M.processInputs[r, p, t, S_v]
           for S_o in M.ProcessOutputsByInput[r, p, t, S_v, S_i]
@@ -1745,7 +1754,7 @@ def MinActivity_Constraint(M, r, p, t):
     except:
       activity_rpt = sum(
           M.V_FlowOutAnnual[r, p, S_i, t, S_v, S_o]
-          for r in reg if '-' not in r
+          for r in reg
           for S_v in M.processVintages[r, p, t]
           for S_i in M.processInputs[r, p, t, S_v]
           for S_o in M.ProcessOutputsByInput[r, p, t, S_v, S_i]