diff --git a/dexomInit.m b/dexomInit.m
index e4b685f..8d79c03 100644
--- a/dexomInit.m
+++ b/dexomInit.m
@@ -28,7 +28,7 @@ function dexomInit(cobraToolboxInitMode, verbose)
     fprintf('This version was tested with Matlab 2015b (CPLEX v12.8), 2018a (CPLEX v12.9),\n');
     fprintf('2018b (CPLEX v12.8) and COBRA Toolbox v3.0.6 on Windows 10.\n');
     
-    paperLink = 'https://doi.org/10.1101/2020.07.17.208918';
+    paperLink = 'https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1008730';
     if usejava('desktop')
         paperLink = ['<a href=\"', paperLink, '\">', paperLink, '</a>'];
     end
diff --git a/src/methods/dexom/dexomCreateMILP.m b/src/methods/dexom/dexomCreateMILP.m
index 051a193..219bb07 100644
--- a/src/methods/dexom/dexomCreateMILP.m
+++ b/src/methods/dexom/dexomCreateMILP.m
@@ -65,6 +65,15 @@
     lb = model.lb;
     ub = model.ub;
 
+    if isempty(RHindex)
+        senseRH = [];
+        numRH = [];
+    end
+    if isempty(RLindex)
+        senseRL = [];
+        numRL = [];
+    end
+
     % Create the A matrix (constraint equations x variables) Ax <=> b
     %
     % Same A matrix as in the original iMAT code +2 rows, one for the numRL
@@ -103,8 +112,9 @@
         A(lastRow+1,i+size(S,2)) = 1;
     end
 
+    lastRow = size(A, 1);
     for i = 1:length(RLindex)
-        A(lastRow+2, i+size(S,2)+length(RHindex)) = 1;
+        A(lastRow+1, i+size(S,2)+length(RHindex)) = 1;
     end
     
     % Alternatively, use fitScore to put a min constraint on the RH+RL
@@ -151,7 +161,7 @@
     csenseFitScore = 'G';
     csense6(1:size(options.exclude,1)) = 'L'; % All less or eq for options.exclude
     csense = [csense1 csense2 csense3 csense4 csense5 senseRH senseRL csenseFitScore csense6];
-
+    
 
     % Creating lb and ub
     lb_y = zeros(2*length(RHindex)+length(RLindex),1);
diff --git a/src/methods/dexomSolveCobraMILP.m b/src/methods/dexomSolveCobraMILP.m
index 584b958..dd28112 100644
--- a/src/methods/dexomSolveCobraMILP.m
+++ b/src/methods/dexomSolveCobraMILP.m
@@ -2,7 +2,9 @@
     global CBT_MILP_SOLVER;
     % Add solver-specific configuration
     if options.useRandomSeed, options.rseed = randi(intmax('int16')); end
-    options.solver = struct;
+    if ~isfield(options,'solver')
+        options.solver = struct;
+    end
     if strcmp(CBT_MILP_SOLVER, 'ibm_cplex')
         options.solver.output.clonelog = -1;
         options.solver.workdir = tempdir();