diff --git a/xbanalysis/problems/gba_problem.py b/xbanalysis/problems/gba_problem.py
index acafcc83d3cf21881d0387d5431a9b0b386fd52a..dee8097ba318e7c0c435f57451aa4f4390fb530c 100644
--- a/xbanalysis/problems/gba_problem.py
+++ b/xbanalysis/problems/gba_problem.py
@@ -26,6 +26,7 @@ def cache_last(func):
 
 # TODO: call get_reaction_rates from solver ?
 # TODO: split get_reaction_rates, get_reaction_jacs, get_reaction_hesses
+# TODO: variable transformation x = exp(y) and solve for unbounded y. start with vector y = log(x)
 
 
 class GbaProblem:
@@ -40,7 +41,9 @@ class GbaProblem:
         :type xba_model: XbaModel
         """
         self.xba_model = xba_model
-        # self.cache = True
+
+        # check to scale density constraint to improve convergence
+        self.scale_density = 1e-4
 
         # metabolites and enzymes in the model
         metab_sids = [s.id for s in xba_model.species.values() if s.sboterm.is_in_branch('SBO:0000247')]
@@ -713,7 +716,7 @@ class GbaProblem:
             density = self.model_params['mws'][self.var_mask_enz].dot(ce)
         else:
             density = self.model_params['mws'].dot(x)
-        return np.array([self.model_params['rho'] - density])
+        return np.array([self.model_params['rho'] - density]) * self.scale_density
 
     @cache_last
     def heq_density_jac(self, x):
@@ -733,7 +736,7 @@ class GbaProblem:
             density_grad[self.var_mask_enz] = -self.model_params['mws'][self.var_mask_enz]
         else:
             density_grad = -self.model_params['mws']
-        return density_grad
+        return density_grad * self.scale_density
 
     # noinspection PyUnusedLocal
     @cache_last
@@ -746,4 +749,4 @@ class GbaProblem:
         :rtype: numpy.ndarray (3D), shape [1,x,x]
         """
         hess = np.zeros((1, self.n_vars, self.n_vars))
-        return np.tril(hess)
+        return np.tril(hess) * self.scale_density