diff --git a/doc/whatsnew.rst b/doc/whatsnew.rst
index d8fb9c9d..59fd9ac7 100644
--- a/doc/whatsnew.rst
+++ b/doc/whatsnew.rst
@@ -26,6 +26,7 @@ New features:
 
 Bug fixes/enhancements:
 
+- fix bug when setting ``param.vary=True`` for a constrained parameter (Issue #859; PR #860)
 - fix bug in reported uncertainties for constrained parameters by better propating uncertainties (Issue #855; PR #856)
 - Coercing of user input data and independent data for ``Model`` to float64 ndarrays is somewhat less aggressive and
   will not increase the precision of numpy ndarrays (see :ref:`model_data_coercion_section` for details). The resulting
diff --git a/examples/doc_builtinmodels_peakmodels.py b/examples/doc_builtinmodels_peakmodels.py
index 7e4a935d..8aa7daa0 100644
--- a/examples/doc_builtinmodels_peakmodels.py
+++ b/examples/doc_builtinmodels_peakmodels.py
@@ -49,9 +49,11 @@
 axes[0].plot(x, out.best_fit, '-', label='Voigt Model\ngamma constrained')
 axes[0].legend()
 
-# free gamma parameter
-pars['gamma'].set(value=0.7, vary=True, expr='')
+# allow the gamma parameter to vary in the fit
+pars['gamma'].vary = True
 out_gamma = mod.fit(y, pars, x=x)
+print(out.fit_report(correl_mode='table'))
+
 axes[1].plot(x, y, '-')
 axes[1].plot(x, out_gamma.best_fit, '-', label='Voigt Model\ngamma unconstrained')
 axes[1].legend()
diff --git a/lmfit/parameter.py b/lmfit/parameter.py
index 224ac55e..d5959979 100644
--- a/lmfit/parameter.py
+++ b/lmfit/parameter.py
@@ -639,7 +639,7 @@ def __init__(self, name, value=None, vary=True, min=-inf, max=inf,
         self.min = min
         self.max = max
         self.brute_step = brute_step
-        self.vary = vary
+        self._vary = vary
         self._expr = expr
         self._expr_ast = None
         self._expr_eval = None
@@ -702,7 +702,7 @@ def set(self, value=None, vary=None, min=None, max=None, expr=None,
 
         """
         if vary is not None:
-            self.vary = vary
+            self._vary = vary
             if vary:
                 self.__set_expression('')
 
@@ -751,13 +751,13 @@ def _init_bounds(self):
 
     def __getstate__(self):
         """Get state for pickle."""
-        return (self.name, self.value, self.vary, self.expr, self.min,
+        return (self.name, self.value, self._vary, self.expr, self.min,
                 self.max, self.brute_step, self.stderr, self.correl,
                 self.init_value, self.user_data)
 
     def __setstate__(self, state):
         """Set state for pickle."""
-        (self.name, _value, self.vary, self.expr, self.min, self.max,
+        (self.name, _value, self._vary, self.expr, self.min, self.max,
          self.brute_step, self.stderr, self.correl, self.init_value,
          self.user_data) = state
         self._expr_ast = None
@@ -772,7 +772,7 @@ def __repr__(self):
         """Return printable representation of a Parameter object."""
         s = []
         sval = f"value={repr(self._getval())}"
-        if not self.vary and self._expr is None:
+        if not self._vary and self._expr is None:
             sval += " (fixed)"
         elif self.stderr is not None:
             sval += f" +/- {self.stderr:.3g}"
@@ -882,6 +882,18 @@ def value(self, val):
         if self._expr_eval is not None:
             self._expr_eval.symtable[self.name] = self._val
 
+    @property
+    def vary(self):
+        """Return whether the parameter is variable"""
+        return self._vary
+
+    @vary.setter
+    def vary(self, val):
+        """Set whether a parameter is varied"""
+        self._vary = val
+        if val:
+            self.__set_expression('')
+
     @property
     def expr(self):
         """Return the mathematical expression used to constrain the value in fit."""
@@ -901,7 +913,7 @@ def __set_expression(self, val):
             val = None
         self._expr = val
         if val is not None:
-            self.vary = False
+            self._vary = False
         if not hasattr(self, '_expr_eval'):
             self._expr_eval = None
         if val is None:
diff --git a/tests/test_parameters.py b/tests/test_parameters.py
index ae39f572..6fea71d8 100644
--- a/tests/test_parameters.py
+++ b/tests/test_parameters.py
@@ -1,7 +1,7 @@
 """Tests for the Parameters class."""
 
-
 from copy import copy, deepcopy
+import os
 import pickle
 
 import numpy as np
@@ -9,6 +9,7 @@
 import pytest
 
 import lmfit
+from lmfit.models import VoigtModel
 
 
 @pytest.fixture
@@ -605,3 +606,31 @@ def test_create_params():
     assert pars1['d'].value == 11
     assert pars1['e'].value == 10000
     assert pars1['e'].brute_step == 4
+
+
+def test_unset_constrained_param():
+    """test 'unsetting' a constrained parameter by
+    just setting `param.vary = True`
+
+    """
+    data = np.loadtxt(os.path.join(os.path.dirname(__file__), '..',
+                                   'examples', 'test_peak.dat'))
+    x = data[:, 0]
+    y = data[:, 1]
+
+    # initial fit
+    mod = VoigtModel()
+    params = mod.guess(y, x=x)
+    out1 = mod.fit(y, params, x=x)
+
+    assert out1.nvarys == 3
+    assert out1.chisqr < 20.0
+
+    # now just gamma to vary
+    params['gamma'].vary = True
+    out2 = mod.fit(y, params, x=x)
+
+    assert out2.nvarys == 4
+    assert out2.chisqr < out1.chisqr
+    assert out2.rsquared > out1.rsquared
+    assert out2.params['gamma'].correl['sigma'] < -0.6