diff --git a/docs/src/api/project.rst b/docs/src/api/project.rst
new file mode 100644
index 00000000..277157b0
--- /dev/null
+++ b/docs/src/api/project.rst
@@ -0,0 +1,7 @@
+Project
+=======
+
+.. autoclass:: easyreflectometry.project.Project
+    :members:
+    :undoc-members:
+    :show-inheritance:
diff --git a/docs/src/tutorials/basic/assemblies_library.rst b/docs/src/tutorials/basic/assemblies_library.rst
index 7e79844f..d72201a9 100644
--- a/docs/src/tutorials/basic/assemblies_library.rst
+++ b/docs/src/tutorials/basic/assemblies_library.rst
@@ -153,8 +153,97 @@ Furthermore, as shown in the `surfactant monolayer tutorial`_ the conformal roug
 
 The use of the :py:class:`SurfactantLayer` in multiple contrast data analysis is shown in a `multiple contrast tutorial`_. 
 
+:py:class:`Bilayer`
+-------------------
+
+The :py:class:`Bilayer` assembly type represents a phospholipid bilayer at an interface.
+It consists of two surfactant layers where one is inverted, creating the structure:
+
+.. code-block:: text
+
+    Head₁ - Tail₁ - Tail₂ - Head₂
+
+This assembly is particularly useful for studying supported lipid bilayers and membrane systems.
+The bilayer comes pre-populated with physically meaningful constraints:
+
+- Both tail layers share the same structural parameters (thickness, area per molecule)
+- Head layers share thickness and area per molecule (different hydration/solvent fraction allowed)
+- A single roughness parameter applies to all interfaces (conformal roughness)
+
+These default constraints can be enabled or disabled as needed for specific analyses.
+
+The creation of a :py:class:`Bilayer` object is shown below.
+
+.. code-block:: python
+
+    from easyreflectometry.sample import Bilayer
+    from easyreflectometry.sample import LayerAreaPerMolecule
+    from easyreflectometry.sample import Material
+
+    # Create materials for solvents
+    d2o = Material(sld=6.36, isld=0.0, name='D2O')
+    air = Material(sld=0.0, isld=0.0, name='Air')
+    
+    # Create head layer (used for front, back head will be auto-created with constraints)
+    head = LayerAreaPerMolecule(
+        molecular_formula='C10H18NO8P',
+        thickness=10.0,
+        solvent=d2o,
+        solvent_fraction=0.3,
+        area_per_molecule=48.2,
+        roughness=3.0,
+        name='DPPC Head'
+    )
+    
+    # Create tail layer (both tail positions will share these parameters)
+    tail = LayerAreaPerMolecule(
+        molecular_formula='C32D64',
+        thickness=16.0,
+        solvent=air,
+        solvent_fraction=0.0,
+        area_per_molecule=48.2,
+        roughness=3.0,
+        name='DPPC Tail'
+    )
+    
+    # Create bilayer with default constraints
+    bilayer = Bilayer(
+        front_head_layer=head,
+        tail_layer=tail,
+        constrain_heads=True,
+        conformal_roughness=True,
+        name='DPPC Bilayer'
+    )
+
+The head layers can have different solvent fractions (hydration) even when constrained, 
+enabling the modeling of asymmetric bilayers at interfaces where the two sides of the 
+bilayer may have different solvent exposure.
+
+The constraints can be controlled at runtime:
+
+.. code-block:: python
+
+    # Disable head constraints to allow different head layer structures
+    bilayer.constrain_heads = False
+    
+    # Disable conformal roughness to allow different roughness values
+    bilayer.conformal_roughness = False
+
+Individual layers can be accessed via properties:
+
+.. code-block:: python
+
+    # Access the four layers
+    bilayer.front_head_layer  # First head layer
+    bilayer.front_tail_layer  # First tail layer  
+    bilayer.back_tail_layer   # Second tail layer (constrained to front tail)
+    bilayer.back_head_layer   # Second head layer
+
+For more detailed examples including simulation and parameter access, see the `bilayer tutorial`_.
+
 
 .. _`simple fitting tutorial`: ../tutorials/simple_fitting.html
 .. _`tutorial`: ../tutorials/repeating.html
 .. _`surfactant monolayer tutorial`: ../tutorials/monolayer.html
-.. _`multiple contrast tutorial`: ../tutorials/multi_contrast.html
\ No newline at end of file
+.. _`multiple contrast tutorial`: ../tutorials/multi_contrast.html
+.. _`bilayer tutorial`: ../tutorials/simulation/bilayer.html
\ No newline at end of file
diff --git a/docs/src/tutorials/simulation/bilayer.ipynb b/docs/src/tutorials/simulation/bilayer.ipynb
new file mode 100644
index 00000000..9725ab1e
--- /dev/null
+++ b/docs/src/tutorials/simulation/bilayer.ipynb
@@ -0,0 +1,542 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "b5afc8c4",
+   "metadata": {},
+   "source": [
+    "# Simulating a Phospholipid Bilayer\n",
+    "\n",
+    "Phospholipid bilayers are fundamental structures in biological membranes and are commonly studied using neutron and X-ray reflectometry.\n",
+    "In this tutorial, we will explore how to use the `Bilayer` assembly in `easyreflectometry` to model a lipid bilayer structure.\n",
+    "\n",
+    "A bilayer consists of two surfactant layers arranged in an inverted configuration:\n",
+    "\n",
+    "```\n",
+    "Head₁ - Tail₁ - Tail₂ - Head₂\n",
+    "```\n",
+    "\n",
+    "The `Bilayer` assembly comes with pre-configured constraints that represent physically meaningful relationships:\n",
+    "- Both tail layers share the same structural parameters\n",
+    "- Head layers share thickness and area per molecule (but can have different hydration)\n",
+    "- Conformal roughness across all interfaces"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f6021005",
+   "metadata": {},
+   "source": [
+    "## Setup\n",
+    "\n",
+    "First, we import the necessary modules and configure matplotlib for inline plotting."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1bfa00f4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "import easyreflectometry\n",
+    "from easyreflectometry.calculators import CalculatorFactory\n",
+    "from easyreflectometry.sample import Bilayer\n",
+    "from easyreflectometry.sample import LayerAreaPerMolecule\n",
+    "from easyreflectometry.sample import Material\n",
+    "from easyreflectometry.sample import Layer\n",
+    "from easyreflectometry.sample import Multilayer\n",
+    "from easyreflectometry.sample import Sample\n",
+    "from easyreflectometry.model import Model\n",
+    "from easyreflectometry.model import PercentageFwhm"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "be41f6f6",
+   "metadata": {},
+   "source": [
+    "## Creating a Bilayer\n",
+    "\n",
+    "We'll create a DPPC (dipalmitoylphosphatidylcholine) bilayer, a common model phospholipid.\n",
+    "\n",
+    "First, let's define the materials for our solvents - D₂O (heavy water) for the head groups and air for the tail groups."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "76bd9056",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Define solvent materials\n",
+    "d2o = Material(sld=6.36, isld=0, name='D2O')\n",
+    "air = Material(sld=0, isld=0, name='Air')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "34a0da6c",
+   "metadata": {},
+   "source": [
+    "### Creating Layer Components\n",
+    "\n",
+    "Now we create the head and tail layers using `LayerAreaPerMolecule`. This approach allows us to define layers based on their chemical formula and area per molecule, which provides a more physically meaningful parameterization."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "54980c8e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Create a head layer for the bilayer\n",
+    "# The head group of DPPC has formula C10H18NO8P\n",
+    "head_layer = LayerAreaPerMolecule(\n",
+    "    molecular_formula='C10H18NO8P',\n",
+    "    thickness=10.0,\n",
+    "    solvent=d2o,\n",
+    "    solvent_fraction=0.3,  # 30% solvent in head region\n",
+    "    area_per_molecule=48.2,\n",
+    "    roughness=3.0,\n",
+    "    name='DPPC Head'\n",
+    ")\n",
+    "\n",
+    "# Create a tail layer for the bilayer\n",
+    "# The tail group of deuterated DPPC has formula C32D64\n",
+    "tail_layer = LayerAreaPerMolecule(\n",
+    "    molecular_formula='C32D64',\n",
+    "    thickness=16.0,\n",
+    "    solvent=air,\n",
+    "    solvent_fraction=0.0,  # No solvent in the tail region\n",
+    "    area_per_molecule=48.2,\n",
+    "    roughness=3.0,\n",
+    "    name='DPPC Tail'\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ce3083b9",
+   "metadata": {},
+   "source": [
+    "### Creating the Bilayer Assembly\n",
+    "\n",
+    "Now we create the `Bilayer` assembly. The bilayer will automatically:\n",
+    "- Create a second tail layer with parameters constrained to the first\n",
+    "- Create a back head layer with thickness and area per molecule constrained to the front head\n",
+    "- Apply conformal roughness across all layers"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4afdd40f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Create the bilayer with default constraints\n",
+    "bilayer = Bilayer(\n",
+    "    front_head_layer=head_layer,\n",
+    "    tail_layer=tail_layer,\n",
+    "    constrain_heads=True,      # Head layers share thickness and area per molecule\n",
+    "    conformal_roughness=True,  # All layers share the same roughness\n",
+    "    name='DPPC Bilayer'\n",
+    ")\n",
+    "\n",
+    "print(bilayer)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "87e39d39",
+   "metadata": {},
+   "source": [
+    "## Exploring the Bilayer Structure\n",
+    "\n",
+    "Let's examine the layer structure of our bilayer."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cc5ffe80",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# The bilayer has 4 layers\n",
+    "print(f'Number of layers: {len(bilayer.layers)}')\n",
+    "print()\n",
+    "\n",
+    "# Access individual layers\n",
+    "print('Layer structure (from front to back):')\n",
+    "print(f'  1. Front Head: {bilayer.front_head_layer.name}')\n",
+    "print(f'  2. Front Tail: {bilayer.front_tail_layer.name}')\n",
+    "print(f'  3. Back Tail:  {bilayer.back_tail_layer.name}')\n",
+    "print(f'  4. Back Head:  {bilayer.back_head_layer.name}')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d9b60931",
+   "metadata": {},
+   "source": [
+    "### Verifying Constraints\n",
+    "\n",
+    "Let's verify that the constraints are working correctly."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a607ca1d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Check tail layer constraints - both tails should have the same parameters\n",
+    "print('Tail layer parameters:')\n",
+    "print(f'  Front tail thickness: {bilayer.front_tail_layer.thickness.value:.2f} Å')\n",
+    "print(f'  Back tail thickness:  {bilayer.back_tail_layer.thickness.value:.2f} Å')\n",
+    "print(f'  Front tail APM: {bilayer.front_tail_layer.area_per_molecule:.2f} Ų')\n",
+    "print(f'  Back tail APM:  {bilayer.back_tail_layer.area_per_molecule:.2f} Ų')\n",
+    "print()\n",
+    "\n",
+    "# Change front tail thickness - back should follow\n",
+    "bilayer.front_tail_layer.thickness.value = 18.0\n",
+    "print('After changing front tail thickness to 18 Å:')\n",
+    "print(f'  Front tail thickness: {bilayer.front_tail_layer.thickness.value:.2f} Å')\n",
+    "print(f'  Back tail thickness:  {bilayer.back_tail_layer.thickness.value:.2f} Å')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ffb84f80",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Check head layer constraints\n",
+    "print('Head layer parameters:')\n",
+    "print(f'  Front head thickness: {bilayer.front_head_layer.thickness.value:.2f} Å')\n",
+    "print(f'  Back head thickness:  {bilayer.back_head_layer.thickness.value:.2f} Å')\n",
+    "print()\n",
+    "\n",
+    "# Hydration (solvent fraction) remains independent\n",
+    "print('Head layer hydration (independent):')\n",
+    "print(f'  Front head solvent fraction: {bilayer.front_head_layer.solvent_fraction:.2f}')\n",
+    "print(f'  Back head solvent fraction:  {bilayer.back_head_layer.solvent_fraction:.2f}')\n",
+    "\n",
+    "# Change back head hydration - front should NOT change\n",
+    "bilayer.back_head_layer.solvent_fraction = 0.5\n",
+    "print()\n",
+    "print('After changing back head solvent fraction to 0.5:')\n",
+    "print(f'  Front head solvent fraction: {bilayer.front_head_layer.solvent_fraction:.2f}')\n",
+    "print(f'  Back head solvent fraction:  {bilayer.back_head_layer.solvent_fraction:.2f}')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "39802839",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Check conformal roughness\n",
+    "print('Roughness values (conformal):')\n",
+    "print(f'  Front head roughness: {bilayer.front_head_layer.roughness.value:.2f} Å')\n",
+    "print(f'  Front tail roughness: {bilayer.front_tail_layer.roughness.value:.2f} Å')\n",
+    "print(f'  Back tail roughness:  {bilayer.back_tail_layer.roughness.value:.2f} Å')\n",
+    "print(f'  Back head roughness:  {bilayer.back_head_layer.roughness.value:.2f} Å')\n",
+    "print()\n",
+    "\n",
+    "# Change front head roughness - all should follow\n",
+    "bilayer.front_head_layer.roughness.value = 4.0\n",
+    "print('After changing front head roughness to 4.0 Å:')\n",
+    "print(f'  Front head roughness: {bilayer.front_head_layer.roughness.value:.2f} Å')\n",
+    "print(f'  Front tail roughness: {bilayer.front_tail_layer.roughness.value:.2f} Å')\n",
+    "print(f'  Back tail roughness:  {bilayer.back_tail_layer.roughness.value:.2f} Å')\n",
+    "print(f'  Back head roughness:  {bilayer.back_head_layer.roughness.value:.2f} Å')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "84f1206b",
+   "metadata": {},
+   "source": [
+    "## Building a Complete Sample\n",
+    "\n",
+    "To simulate reflectometry, we need to create a complete sample with sub- and super-phases."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1e509350",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Reset bilayer parameters\n",
+    "bilayer.front_head_layer.roughness.value = 3.0\n",
+    "bilayer.front_tail_layer.thickness.value = 16.0\n",
+    "bilayer.back_head_layer.solvent_fraction = 0.3\n",
+    "\n",
+    "# Create superphase (air) and subphase (silicon substrate with oxide layer)\n",
+    "vacuum = Material(sld=0, isld=0, name='Vacuum')\n",
+    "superphase = Layer(material=vacuum, thickness=0, roughness=0, name='Vacuum Superphase')\n",
+    "\n",
+    "si = Material(sld=2.047, isld=0, name='Si')\n",
+    "sio2 = Material(sld=3.47, isld=0, name='SiO2')\n",
+    "si_layer = Layer(material=si, thickness=0, roughness=0, name='Si')\n",
+    "sio2_layer = Layer(material=sio2, thickness=15, roughness=3, name='SiO2')\n",
+    "\n",
+    "# D2O subphase\n",
+    "d2o_layer = Layer(material=d2o, thickness=0, roughness=3, name='D2O Subphase')\n",
+    "\n",
+    "# Create sample: superphase -> bilayer -> SiO2 -> Si\n",
+    "sample = Sample(\n",
+    "    Multilayer(superphase, name='Superphase'),\n",
+    "    bilayer,\n",
+    "    Multilayer(d2o_layer, name='D2O'),\n",
+    "    Multilayer([sio2_layer, si_layer], name='Substrate'),\n",
+    "    name='Bilayer on Si/SiO2'\n",
+    ")\n",
+    "\n",
+    "print(sample)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "0ce27fe5",
+   "metadata": {},
+   "source": [
+    "## Simulating Reflectivity\n",
+    "\n",
+    "Now we can simulate the reflectometry profile for our bilayer sample."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "58538a77",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Create the model\n",
+    "model = Model(\n",
+    "    sample=sample,\n",
+    "    scale=1.0,\n",
+    "    background=1e-7,\n",
+    "    resolution_function=PercentageFwhm(5),\n",
+    "    name='Bilayer Model'\n",
+    ")\n",
+    "\n",
+    "# Set up the calculator\n",
+    "interface = CalculatorFactory()\n",
+    "model.interface = interface\n",
+    "\n",
+    "# Generate Q values\n",
+    "q = np.linspace(0.005, 0.3, 500)\n",
+    "\n",
+    "# Calculate reflectometry\n",
+    "reflectivity = model.interface().reflectity_profile(q, model.unique_name)\n",
+    "\n",
+    "# Plot\n",
+    "plt.figure(figsize=(10, 6))\n",
+    "plt.semilogy(q, reflectivity, 'b-', linewidth=2, label='Bilayer')\n",
+    "plt.xlabel('Q (Å⁻¹)')\n",
+    "plt.ylabel('Reflectivity')\n",
+    "plt.title('Simulated Reflectometry of DPPC Bilayer')\n",
+    "plt.legend()\n",
+    "plt.grid(True, alpha=0.3)\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "06723f8f",
+   "metadata": {},
+   "source": [
+    "## Scattering Length Density Profile\n",
+    "\n",
+    "Let's also visualize the SLD profile of our bilayer structure."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a61981e0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Get SLD profile\n",
+    "z, sld = model.interface().sld_profile(model.unique_name)\n",
+    "\n",
+    "plt.figure(figsize=(10, 6))\n",
+    "plt.plot(z, sld, 'b-', linewidth=2)\n",
+    "plt.xlabel('Distance from interface (Å)')\n",
+    "plt.ylabel('SLD (10⁻⁶ Å⁻²)')\n",
+    "plt.title('SLD Profile of DPPC Bilayer on Si/SiO2')\n",
+    "plt.grid(True, alpha=0.3)\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "504f5fbe",
+   "metadata": {},
+   "source": [
+    "## Modifying Constraints\n",
+    "\n",
+    "The bilayer constraints can be modified at runtime. Let's see how disabling conformal roughness affects the model."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "12bffbbb",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Disable conformal roughness\n",
+    "bilayer.conformal_roughness = False\n",
+    "\n",
+    "# Now we can set different roughness values for each layer\n",
+    "bilayer.front_head_layer.roughness.value = 2.0\n",
+    "bilayer.front_tail_layer.roughness.value = 1.5\n",
+    "bilayer.back_tail_layer.roughness.value = 1.5\n",
+    "bilayer.back_head_layer.roughness.value = 4.0\n",
+    "\n",
+    "print('Individual roughness values after disabling conformal roughness:')\n",
+    "print(f'  Front head: {bilayer.front_head_layer.roughness.value:.2f} Å')\n",
+    "print(f'  Front tail: {bilayer.front_tail_layer.roughness.value:.2f} Å')\n",
+    "print(f'  Back tail:  {bilayer.back_tail_layer.roughness.value:.2f} Å')\n",
+    "print(f'  Back head:  {bilayer.back_head_layer.roughness.value:.2f} Å')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5fd26d95",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Compare reflectometry with different roughness configurations\n",
+    "reflectivity_variable_roughness = model.interface().reflectity_profile(q, model.unique_name)\n",
+    "\n",
+    "# Reset to conformal roughness\n",
+    "bilayer.conformal_roughness = True\n",
+    "bilayer.front_head_layer.roughness.value = 3.0\n",
+    "reflectivity_conformal = model.interface().reflectity_profile(q, model.unique_name)\n",
+    "\n",
+    "plt.figure(figsize=(10, 6))\n",
+    "plt.semilogy(q, reflectivity_conformal, 'b-', linewidth=2, label='Conformal roughness (3.0 Å)')\n",
+    "plt.semilogy(q, reflectivity_variable_roughness, 'r--', linewidth=2, label='Variable roughness')\n",
+    "plt.xlabel('Q (Å⁻¹)')\n",
+    "plt.ylabel('Reflectivity')\n",
+    "plt.title('Effect of Roughness Configuration on Reflectometry')\n",
+    "plt.legend()\n",
+    "plt.grid(True, alpha=0.3)\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "955e1320",
+   "metadata": {},
+   "source": [
+    "## Asymmetric Hydration\n",
+    "\n",
+    "One of the key features of the `Bilayer` assembly is the ability to model asymmetric hydration - where the two sides of the bilayer have different solvent exposure."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "993ad6d4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Create a bilayer with asymmetric hydration\n",
+    "# This is common in supported bilayers where one side is in contact with a substrate\n",
+    "\n",
+    "# Low hydration on front (substrate side)\n",
+    "bilayer.front_head_layer.solvent_fraction = 0.1\n",
+    "\n",
+    "# Higher hydration on back (solution side)\n",
+    "bilayer.back_head_layer.solvent_fraction = 0.4\n",
+    "\n",
+    "print('Asymmetric hydration:')\n",
+    "print(f'  Front head solvent fraction: {bilayer.front_head_layer.solvent_fraction:.2f}')\n",
+    "print(f'  Back head solvent fraction:  {bilayer.back_head_layer.solvent_fraction:.2f}')\n",
+    "\n",
+    "# Calculate new reflectometry\n",
+    "reflectivity_asymmetric = model.interface().reflectity_profile(q, model.unique_name)\n",
+    "\n",
+    "# Reset to symmetric for comparison\n",
+    "bilayer.front_head_layer.solvent_fraction = 0.3\n",
+    "bilayer.back_head_layer.solvent_fraction = 0.3\n",
+    "reflectivity_symmetric = model.interface().reflectity_profile(q, model.unique_name)\n",
+    "\n",
+    "plt.figure(figsize=(10, 6))\n",
+    "plt.semilogy(q, reflectivity_symmetric, 'b-', linewidth=2, label='Symmetric hydration (0.3/0.3)')\n",
+    "plt.semilogy(q, reflectivity_asymmetric, 'r--', linewidth=2, label='Asymmetric hydration (0.1/0.4)')\n",
+    "plt.xlabel('Q (Å⁻¹)')\n",
+    "plt.ylabel('Reflectivity')\n",
+    "plt.title('Effect of Asymmetric Hydration on Reflectometry')\n",
+    "plt.legend()\n",
+    "plt.grid(True, alpha=0.3)\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "df3d351f",
+   "metadata": {},
+   "source": [
+    "## Summary\n",
+    "\n",
+    "In this tutorial, we have explored the `Bilayer` assembly in `easyreflectometry`:\n",
+    "\n",
+    "1. **Creating a bilayer**: Using `LayerAreaPerMolecule` components for head and tail layers\n",
+    "2. **Understanding constraints**: \n",
+    "   - Tail layers share all structural parameters\n",
+    "   - Head layers share thickness and area per molecule (hydration is independent)\n",
+    "   - Conformal roughness applies to all layers by default\n",
+    "3. **Building a sample**: Combining the bilayer with sub- and super-phases\n",
+    "4. **Simulating reflectometry**: Using the calculator interface to generate reflectivity profiles\n",
+    "5. **Modifying constraints**: Disabling conformal roughness for more complex models\n",
+    "6. **Asymmetric hydration**: Modeling supported bilayers with different solvent exposure on each side\n",
+    "\n",
+    "The `Bilayer` assembly provides a convenient way to model phospholipid bilayers with physically meaningful constraints, reducing the number of free parameters while maintaining flexibility for complex systems."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "era",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.11"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/docs/src/tutorials/simulation/simulation.rst b/docs/src/tutorials/simulation/simulation.rst
index 4f187b05..5cf9ced4 100644
--- a/docs/src/tutorials/simulation/simulation.rst
+++ b/docs/src/tutorials/simulation/simulation.rst
@@ -6,5 +6,6 @@ These are basic simulation examples using the :py:mod:`easyreflectometry` librar
 .. toctree:: 
    :maxdepth: 1
    
+   bilayer.ipynb
    magnetism.ipynb
    resolution_functions.ipynb
\ No newline at end of file
diff --git a/src/easyreflectometry/sample/__init__.py b/src/easyreflectometry/sample/__init__.py
index 2012cd44..4991b975 100644
--- a/src/easyreflectometry/sample/__init__.py
+++ b/src/easyreflectometry/sample/__init__.py
@@ -1,4 +1,5 @@
 from .assemblies.base_assembly import BaseAssembly
+from .assemblies.bilayer import Bilayer
 from .assemblies.gradient_layer import GradientLayer
 from .assemblies.multilayer import Multilayer
 from .assemblies.repeating_multilayer import RepeatingMultilayer
@@ -15,6 +16,7 @@
 
 __all__ = (
     'BaseAssembly',
+    'Bilayer',
     'GradientLayer',
     'Layer',
     'LayerAreaPerMolecule',
diff --git a/src/easyreflectometry/sample/assemblies/bilayer.py b/src/easyreflectometry/sample/assemblies/bilayer.py
new file mode 100644
index 00000000..888c4d1d
--- /dev/null
+++ b/src/easyreflectometry/sample/assemblies/bilayer.py
@@ -0,0 +1,425 @@
+from __future__ import annotations
+
+from typing import Optional
+
+from easyscience import global_object
+from easyscience.variable import Parameter
+
+from ..collections.layer_collection import LayerCollection
+from ..elements.layers.layer_area_per_molecule import LayerAreaPerMolecule
+from ..elements.materials.material import Material
+from .base_assembly import BaseAssembly
+
+
+class Bilayer(BaseAssembly):
+    """A lipid bilayer consisting of two surfactant layers where one is inverted.
+
+    The bilayer structure is: Front Head - Front Tail - Back Tail - Back Head
+
+    This assembly comes pre-populated with physically meaningful constraints:
+    - Both tail layers share the same structural parameters (thickness, area per molecule)
+    - Head layers share thickness and area per molecule (different hydration/solvent_fraction allowed)
+    - A single roughness parameter applies to all interfaces (conformal roughness)
+
+    More information about the usage of this assembly is available in the
+    `bilayer documentation`_
+
+    .. _`bilayer documentation`: ../sample/assemblies_library.html#bilayer
+    """
+
+    def __init__(
+        self,
+        front_head_layer: Optional[LayerAreaPerMolecule] = None,
+        tail_layer: Optional[LayerAreaPerMolecule] = None,
+        back_head_layer: Optional[LayerAreaPerMolecule] = None,
+        name: str = 'EasyBilayer',
+        unique_name: Optional[str] = None,
+        constrain_heads: bool = True,
+        conformal_roughness: bool = True,
+        interface=None,
+    ):
+        """Constructor.
+
+        :param front_head_layer: Layer representing the front head part of the bilayer.
+        :param tail_layer: Layer representing the tail part of the bilayer. A second tail
+            layer is created internally with parameters constrained to this one.
+        :param back_head_layer: Layer representing the back head part of the bilayer.
+        :param name: Name for bilayer, defaults to 'EasyBilayer'.
+        :param constrain_heads: Constrain head layer thickness and area per molecule, defaults to `True`.
+        :param conformal_roughness: Constrain roughness to be the same for all layers, defaults to `True`.
+        :param interface: Calculator interface, defaults to `None`.
+        """
+        # Generate unique name for nested objects
+        if unique_name is None:
+            unique_name = global_object.generate_unique_name(self.__class__.__name__)
+
+        # Create default front head layer if not provided
+        if front_head_layer is None:
+            d2o_front = Material(
+                sld=6.36,
+                isld=0,
+                name='D2O',
+                unique_name=unique_name + '_MaterialFrontHead',
+                interface=interface,
+            )
+            front_head_layer = LayerAreaPerMolecule(
+                molecular_formula='C10H18NO8P',
+                thickness=10.0,
+                solvent=d2o_front,
+                solvent_fraction=0.2,
+                area_per_molecule=48.2,
+                roughness=3.0,
+                name='DPPC Head Front',
+                unique_name=unique_name + '_LayerAreaPerMoleculeFrontHead',
+                interface=interface,
+            )
+
+        # Create default tail layer if not provided
+        if tail_layer is None:
+            air = Material(
+                sld=0,
+                isld=0,
+                name='Air',
+                unique_name=unique_name + '_MaterialTail',
+                interface=interface,
+            )
+            tail_layer = LayerAreaPerMolecule(
+                molecular_formula='C32D64',
+                thickness=16.0,
+                solvent=air,
+                solvent_fraction=0.0,
+                area_per_molecule=48.2,
+                roughness=3.0,
+                name='DPPC Tail',
+                unique_name=unique_name + '_LayerAreaPerMoleculeTail',
+                interface=interface,
+            )
+
+        # Create second tail layer with same parameters as the first
+        # This will be constrained to the first tail layer
+        air_back = Material(
+            sld=0,
+            isld=0,
+            name='Air',
+            unique_name=unique_name + '_MaterialBackTail',
+            interface=interface,
+        )
+        back_tail_layer = LayerAreaPerMolecule(
+            molecular_formula=tail_layer.molecular_formula,
+            thickness=tail_layer.thickness.value,
+            solvent=air_back,
+            solvent_fraction=tail_layer.solvent_fraction,
+            area_per_molecule=tail_layer.area_per_molecule,
+            roughness=tail_layer.roughness.value,
+            name=tail_layer.name + ' Back',
+            unique_name=unique_name + '_LayerAreaPerMoleculeBackTail',
+            interface=interface,
+        )
+
+        # Create default back head layer if not provided
+        if back_head_layer is None:
+            d2o_back = Material(
+                sld=6.36,
+                isld=0,
+                name='D2O',
+                unique_name=unique_name + '_MaterialBackHead',
+                interface=interface,
+            )
+            back_head_layer = LayerAreaPerMolecule(
+                molecular_formula='C10H18NO8P',
+                thickness=10.0,
+                solvent=d2o_back,
+                solvent_fraction=0.2,
+                area_per_molecule=48.2,
+                roughness=3.0,
+                name='DPPC Head Back',
+                unique_name=unique_name + '_LayerAreaPerMoleculeBackHead',
+                interface=interface,
+            )
+
+        # Store the front tail layer reference for constraint setup
+        self._front_tail_layer = tail_layer
+        self._back_tail_layer = back_tail_layer
+
+        # Create layer collection: front_head, front_tail, back_tail, back_head
+        bilayer_layers = LayerCollection(
+            front_head_layer,
+            tail_layer,
+            back_tail_layer,
+            back_head_layer,
+            name='Layers',
+            unique_name=unique_name + '_LayerCollection',
+            interface=interface,
+        )
+
+        super().__init__(
+            name=name,
+            unique_name=unique_name,
+            type='Bilayer',
+            layers=bilayer_layers,
+            interface=interface,
+        )
+
+        self.interface = interface
+        self._conformal_roughness = False
+        self._constrain_heads = False
+        self._tail_constraints_setup = False
+
+        # Setup tail layer constraints (back tail depends on front tail)
+        self._setup_tail_constraints()
+
+        # Apply head constraints if requested
+        if constrain_heads:
+            self.constrain_heads = True
+
+        # Apply conformal roughness if requested
+        if conformal_roughness:
+            self.conformal_roughness = True
+
+    def _setup_tail_constraints(self) -> None:
+        """Setup constraints so back tail layer parameters depend on front tail layer."""
+        front_tail = self._front_tail_layer
+        back_tail = self._back_tail_layer
+
+        # Constrain thickness
+        back_tail.thickness.make_dependent_on(
+            dependency_expression='a',
+            dependency_map={'a': front_tail.thickness},
+        )
+
+        # Constrain area per molecule
+        back_tail._area_per_molecule.make_dependent_on(
+            dependency_expression='a',
+            dependency_map={'a': front_tail._area_per_molecule},
+        )
+
+        # Constrain solvent fraction
+        back_tail.material._fraction.make_dependent_on(
+            dependency_expression='a',
+            dependency_map={'a': front_tail.material._fraction},
+        )
+
+        self._tail_constraints_setup = True
+
+    @property
+    def front_head_layer(self) -> Optional[LayerAreaPerMolecule]:
+        """Get the front head layer of the bilayer."""
+        return self.layers[0]
+
+    @front_head_layer.setter
+    def front_head_layer(self, layer: LayerAreaPerMolecule) -> None:
+        """Set the front head layer of the bilayer."""
+        self.layers[0] = layer
+
+    @property
+    def front_tail_layer(self) -> Optional[LayerAreaPerMolecule]:
+        """Get the front tail layer of the bilayer."""
+        return self.layers[1]
+
+    @property
+    def tail_layer(self) -> Optional[LayerAreaPerMolecule]:
+        """Get the tail layer (alias for front_tail_layer for serialization compatibility)."""
+        return self.front_tail_layer
+
+    @property
+    def back_tail_layer(self) -> Optional[LayerAreaPerMolecule]:
+        """Get the back tail layer of the bilayer."""
+        return self.layers[2]
+
+    @property
+    def back_head_layer(self) -> Optional[LayerAreaPerMolecule]:
+        """Get the back head layer of the bilayer."""
+        return self.layers[3]
+
+    @back_head_layer.setter
+    def back_head_layer(self, layer: LayerAreaPerMolecule) -> None:
+        """Set the back head layer of the bilayer."""
+        self.layers[3] = layer
+
+    @property
+    def constrain_heads(self) -> bool:
+        """Get the head layer constraint status."""
+        return self._constrain_heads
+
+    @constrain_heads.setter
+    def constrain_heads(self, status: bool) -> None:
+        """Set the status for head layer constraints.
+
+        When enabled, the back head layer thickness and area per molecule
+        are constrained to match the front head layer. Solvent fraction
+        (hydration) remains independent.
+
+        :param status: Boolean for the constraint status.
+        """
+        if status:
+            self._enable_head_constraints()
+        else:
+            self._disable_head_constraints()
+        self._constrain_heads = status
+
+    def _enable_head_constraints(self) -> None:
+        """Enable head layer constraints (thickness and area per molecule)."""
+        front_head = self.front_head_layer
+        back_head = self.back_head_layer
+
+        # Constrain thickness
+        back_head.thickness.make_dependent_on(
+            dependency_expression='a',
+            dependency_map={'a': front_head.thickness},
+        )
+
+        # Constrain area per molecule
+        back_head._area_per_molecule.make_dependent_on(
+            dependency_expression='a',
+            dependency_map={'a': front_head._area_per_molecule},
+        )
+
+    def _disable_head_constraints(self) -> None:
+        """Disable head layer constraints."""
+        self.back_head_layer.thickness.make_independent()
+        self.back_head_layer._area_per_molecule.make_independent()
+
+    @property
+    def conformal_roughness(self) -> bool:
+        """Get the roughness constraint status."""
+        return self._conformal_roughness
+
+    @conformal_roughness.setter
+    def conformal_roughness(self, status: bool) -> None:
+        """Set the status for conformal roughness.
+
+        When enabled, all layers share the same roughness parameter
+        (controlled by the front head layer).
+
+        :param status: Boolean for the constraint status.
+        """
+        if status:
+            self._setup_roughness_constraints()
+            self._enable_roughness_constraints()
+        else:
+            if self._roughness_constraints_setup:
+                self._disable_roughness_constraints()
+        self._conformal_roughness = status
+
+    def constrain_solvent_roughness(self, solvent_roughness: Parameter) -> None:
+        """Add the constraint to the solvent roughness.
+
+        :param solvent_roughness: The solvent roughness parameter.
+        """
+        if not self.conformal_roughness:
+            raise ValueError('Roughness must be conformal to use this function.')
+        solvent_roughness.value = self.front_head_layer.roughness.value
+        solvent_roughness.make_dependent_on(
+            dependency_expression='a',
+            dependency_map={'a': self.front_head_layer.roughness},
+        )
+
+    def constrain_multiple_contrast(
+        self,
+        another_contrast: Bilayer,
+        front_head_thickness: bool = True,
+        back_head_thickness: bool = True,
+        tail_thickness: bool = True,
+        front_head_area_per_molecule: bool = True,
+        back_head_area_per_molecule: bool = True,
+        tail_area_per_molecule: bool = True,
+        front_head_fraction: bool = True,
+        back_head_fraction: bool = True,
+        tail_fraction: bool = True,
+    ) -> None:
+        """Constrain structural parameters between bilayer objects.
+
+        :param another_contrast: The bilayer to constrain to.
+        :param front_head_thickness: Constrain front head thickness.
+        :param back_head_thickness: Constrain back head thickness.
+        :param tail_thickness: Constrain tail thickness.
+        :param front_head_area_per_molecule: Constrain front head area per molecule.
+        :param back_head_area_per_molecule: Constrain back head area per molecule.
+        :param tail_area_per_molecule: Constrain tail area per molecule.
+        :param front_head_fraction: Constrain front head solvent fraction.
+        :param back_head_fraction: Constrain back head solvent fraction.
+        :param tail_fraction: Constrain tail solvent fraction.
+        """
+        if front_head_thickness:
+            self.front_head_layer.thickness.make_dependent_on(
+                dependency_expression='a',
+                dependency_map={'a': another_contrast.front_head_layer.thickness},
+            )
+
+        if back_head_thickness:
+            self.back_head_layer.thickness.make_dependent_on(
+                dependency_expression='a',
+                dependency_map={'a': another_contrast.back_head_layer.thickness},
+            )
+
+        if tail_thickness:
+            self.front_tail_layer.thickness.make_dependent_on(
+                dependency_expression='a',
+                dependency_map={'a': another_contrast.front_tail_layer.thickness},
+            )
+
+        if front_head_area_per_molecule:
+            self.front_head_layer._area_per_molecule.make_dependent_on(
+                dependency_expression='a',
+                dependency_map={'a': another_contrast.front_head_layer._area_per_molecule},
+            )
+
+        if back_head_area_per_molecule:
+            self.back_head_layer._area_per_molecule.make_dependent_on(
+                dependency_expression='a',
+                dependency_map={'a': another_contrast.back_head_layer._area_per_molecule},
+            )
+
+        if tail_area_per_molecule:
+            self.front_tail_layer._area_per_molecule.make_dependent_on(
+                dependency_expression='a',
+                dependency_map={'a': another_contrast.front_tail_layer._area_per_molecule},
+            )
+
+        if front_head_fraction:
+            self.front_head_layer.material._fraction.make_dependent_on(
+                dependency_expression='a',
+                dependency_map={'a': another_contrast.front_head_layer.material._fraction},
+            )
+
+        if back_head_fraction:
+            self.back_head_layer.material._fraction.make_dependent_on(
+                dependency_expression='a',
+                dependency_map={'a': another_contrast.back_head_layer.material._fraction},
+            )
+
+        if tail_fraction:
+            self.front_tail_layer.material._fraction.make_dependent_on(
+                dependency_expression='a',
+                dependency_map={'a': another_contrast.front_tail_layer.material._fraction},
+            )
+
+    @property
+    def _dict_repr(self) -> dict:
+        """A simplified dict representation."""
+        return {
+            self.name: {
+                'front_head_layer': self.front_head_layer._dict_repr,
+                'front_tail_layer': self.front_tail_layer._dict_repr,
+                'back_tail_layer': self.back_tail_layer._dict_repr,
+                'back_head_layer': self.back_head_layer._dict_repr,
+                'constrain_heads': self.constrain_heads,
+                'conformal_roughness': self.conformal_roughness,
+            }
+        }
+
+    def as_dict(self, skip: Optional[list[str]] = None) -> dict:
+        """Produce a cleaned dict using a custom as_dict method.
+
+        The resulting dict matches the parameters in __init__
+
+        :param skip: List of keys to skip, defaults to `None`.
+        """
+        this_dict = super().as_dict(skip=skip)
+        this_dict['front_head_layer'] = self.front_head_layer.as_dict(skip=skip)
+        this_dict['tail_layer'] = self.front_tail_layer.as_dict(skip=skip)
+        this_dict['back_head_layer'] = self.back_head_layer.as_dict(skip=skip)
+        this_dict['constrain_heads'] = self.constrain_heads
+        this_dict['conformal_roughness'] = self.conformal_roughness
+        del this_dict['layers']
+        return this_dict
diff --git a/tests/sample/assemblies/test_bilayer.py b/tests/sample/assemblies/test_bilayer.py
new file mode 100644
index 00000000..a17ae0b0
--- /dev/null
+++ b/tests/sample/assemblies/test_bilayer.py
@@ -0,0 +1,410 @@
+"""
+Tests for Bilayer class module
+"""
+
+__author__ = 'github.com/easyscience'
+__version__ = '0.0.1'
+
+
+from easyscience import global_object
+
+from easyreflectometry.sample.assemblies.bilayer import Bilayer
+from easyreflectometry.sample.elements.layers.layer import Layer
+from easyreflectometry.sample.elements.layers.layer_area_per_molecule import LayerAreaPerMolecule
+from easyreflectometry.sample.elements.materials.material import Material
+
+
+class TestBilayer:
+    def test_default(self):
+        """Test default bilayer creation with expected structure."""
+        p = Bilayer()
+        assert p.name == 'EasyBilayer'
+        assert p._type == 'Bilayer'
+
+        # Check layer count
+        assert len(p.layers) == 4
+
+        # Check layer order: front_head, front_tail, back_tail, back_head
+        assert p.layers[0].name == 'DPPC Head Front'
+        assert p.front_head_layer.name == 'DPPC Head Front'
+
+        assert p.layers[1].name == 'DPPC Tail'
+        assert p.front_tail_layer.name == 'DPPC Tail'
+
+        assert p.layers[2].name == 'DPPC Tail Back'
+        assert p.back_tail_layer.name == 'DPPC Tail Back'
+
+        assert p.layers[3].name == 'DPPC Head Back'
+        assert p.back_head_layer.name == 'DPPC Head Back'
+
+    def test_default_constraints_enabled(self):
+        """Test that default bilayer has constraints enabled."""
+        p = Bilayer()
+
+        # Default should have conformal roughness and head constraints
+        assert p.conformal_roughness is True
+        assert p.constrain_heads is True
+
+    def test_layer_structure(self):
+        """Verify 4 layers in correct order."""
+        p = Bilayer()
+
+        assert p.front_head_layer is p.layers[0]
+        assert p.front_tail_layer is p.layers[1]
+        assert p.back_tail_layer is p.layers[2]
+        assert p.back_head_layer is p.layers[3]
+
+    def test_custom_layers(self):
+        """Test creation with custom head/tail layers."""
+        d2o = Material(sld=6.36, isld=0, name='D2O')
+        air = Material(sld=0, isld=0, name='Air')
+
+        front_head = LayerAreaPerMolecule(
+            molecular_formula='C10H18NO8P',
+            thickness=12.0,
+            solvent=d2o,
+            solvent_fraction=0.3,
+            area_per_molecule=50.0,
+            roughness=2.0,
+            name='Custom Front Head',
+        )
+        tail = LayerAreaPerMolecule(
+            molecular_formula='C32D64',
+            thickness=18.0,
+            solvent=air,
+            solvent_fraction=0.0,
+            area_per_molecule=50.0,
+            roughness=2.0,
+            name='Custom Tail',
+        )
+        back_head = LayerAreaPerMolecule(
+            molecular_formula='C10H18NO8P',
+            thickness=12.0,
+            solvent=d2o,
+            solvent_fraction=0.4,  # Different hydration
+            area_per_molecule=50.0,
+            roughness=2.0,
+            name='Custom Back Head',
+        )
+
+        p = Bilayer(
+            front_head_layer=front_head,
+            tail_layer=tail,
+            back_head_layer=back_head,
+            name='Custom Bilayer',
+        )
+
+        assert p.name == 'Custom Bilayer'
+        assert p.front_head_layer.name == 'Custom Front Head'
+        assert p.front_tail_layer.name == 'Custom Tail'
+        assert p.back_head_layer.name == 'Custom Back Head'
+        assert p.front_head_layer.thickness.value == 12.0
+        assert p.front_tail_layer.thickness.value == 18.0
+
+    def test_tail_layers_linked(self):
+        """Test that both tail layers share parameters."""
+        p = Bilayer()
+
+        # Initial values should match
+        assert p.front_tail_layer.thickness.value == p.back_tail_layer.thickness.value
+        assert p.front_tail_layer.area_per_molecule == p.back_tail_layer.area_per_molecule
+
+        # Change front tail thickness - back tail should follow
+        p.front_tail_layer.thickness.value = 20.0
+        assert p.front_tail_layer.thickness.value == 20.0
+        assert p.back_tail_layer.thickness.value == 20.0
+
+        # Change front tail area per molecule - back tail should follow
+        p.front_tail_layer._area_per_molecule.value = 55.0
+        assert p.front_tail_layer.area_per_molecule == 55.0
+        assert p.back_tail_layer.area_per_molecule == 55.0
+
+    def test_constrain_heads_enabled(self):
+        """Test head thickness/area constraint when enabled."""
+        p = Bilayer(constrain_heads=True)
+
+        # Change front head thickness - back head should follow
+        p.front_head_layer.thickness.value = 15.0
+        assert p.front_head_layer.thickness.value == 15.0
+        assert p.back_head_layer.thickness.value == 15.0
+
+        # Change front head area per molecule - back head should follow
+        p.front_head_layer._area_per_molecule.value = 60.0
+        assert p.front_head_layer.area_per_molecule == 60.0
+        assert p.back_head_layer.area_per_molecule == 60.0
+
+    def test_constrain_heads_disabled(self):
+        """Test heads are independent when constraint disabled."""
+        p = Bilayer(constrain_heads=False)
+
+        # Set different values for front and back heads
+        p.front_head_layer.thickness.value = 15.0
+        p.back_head_layer.thickness.value = 12.0
+
+        assert p.front_head_layer.thickness.value == 15.0
+        assert p.back_head_layer.thickness.value == 12.0
+
+    def test_constrain_heads_toggle(self):
+        """Test enabling/disabling head constraints at runtime."""
+        p = Bilayer(constrain_heads=False)
+
+        # Set different values
+        p.front_head_layer.thickness.value = 15.0
+        p.back_head_layer.thickness.value = 12.0
+
+        # Enable constraint - back head should match front head
+        p.constrain_heads = True
+        assert p.constrain_heads is True
+
+        # Change front head - back should follow
+        p.front_head_layer.thickness.value = 20.0
+        assert p.back_head_layer.thickness.value == 20.0
+
+        # Disable constraint
+        p.constrain_heads = False
+        assert p.constrain_heads is False
+
+        # Now they can be independent
+        p.back_head_layer.thickness.value = 18.0
+        assert p.front_head_layer.thickness.value == 20.0
+        assert p.back_head_layer.thickness.value == 18.0
+
+    def test_head_hydration_independent(self):
+        """Test that head hydrations remain independent even with constraints."""
+        p = Bilayer(constrain_heads=True)
+
+        # Set different solvent fractions
+        p.front_head_layer.solvent_fraction = 0.3
+        p.back_head_layer.solvent_fraction = 0.5
+
+        # They should remain independent
+        assert p.front_head_layer.solvent_fraction == 0.3
+        assert p.back_head_layer.solvent_fraction == 0.5
+
+    def test_conformal_roughness_enabled(self):
+        """Test all roughnesses are linked when conformal roughness enabled."""
+        p = Bilayer(conformal_roughness=True)
+
+        # Change front head roughness - all should follow
+        p.front_head_layer.roughness.value = 5.0
+        assert p.front_head_layer.roughness.value == 5.0
+        assert p.front_tail_layer.roughness.value == 5.0
+        assert p.back_tail_layer.roughness.value == 5.0
+        assert p.back_head_layer.roughness.value == 5.0
+
+    def test_conformal_roughness_disabled(self):
+        """Test roughnesses are independent when conformal roughness disabled."""
+        p = Bilayer(conformal_roughness=False)
+
+        # Set different roughnesses
+        p.front_head_layer.roughness.value = 2.0
+        p.front_tail_layer.roughness.value = 3.0
+        p.back_tail_layer.roughness.value = 4.0
+        p.back_head_layer.roughness.value = 5.0
+
+        assert p.front_head_layer.roughness.value == 2.0
+        assert p.front_tail_layer.roughness.value == 3.0
+        assert p.back_tail_layer.roughness.value == 4.0
+        assert p.back_head_layer.roughness.value == 5.0
+
+    def test_conformal_roughness_toggle(self):
+        """Test enabling/disabling conformal roughness at runtime."""
+        p = Bilayer(conformal_roughness=False)
+
+        # Set different values
+        p.front_head_layer.roughness.value = 2.0
+        p.back_head_layer.roughness.value = 5.0
+
+        # Enable conformal roughness
+        p.conformal_roughness = True
+        assert p.conformal_roughness is True
+
+        # Change front head - all should follow
+        p.front_head_layer.roughness.value = 4.0
+        assert p.front_tail_layer.roughness.value == 4.0
+        assert p.back_tail_layer.roughness.value == 4.0
+        assert p.back_head_layer.roughness.value == 4.0
+
+        # Disable conformal roughness
+        p.conformal_roughness = False
+        assert p.conformal_roughness is False
+
+    def test_get_set_front_head_layer(self):
+        """Test getting and setting front head layer."""
+        p = Bilayer()
+        new_layer = LayerAreaPerMolecule(
+            molecular_formula='C8H16NO6P',
+            thickness=8.0,
+            name='New Front Head',
+        )
+
+        p.front_head_layer = new_layer
+
+        assert p.front_head_layer == new_layer
+        assert p.layers[0] == new_layer
+
+    def test_get_set_back_head_layer(self):
+        """Test getting and setting back head layer."""
+        p = Bilayer()
+        new_layer = LayerAreaPerMolecule(
+            molecular_formula='C8H16NO6P',
+            thickness=8.0,
+            name='New Back Head',
+        )
+
+        p.back_head_layer = new_layer
+
+        assert p.back_head_layer == new_layer
+        assert p.layers[3] == new_layer
+
+    def test_dict_repr(self):
+        """Test dictionary representation."""
+        p = Bilayer()
+
+        dict_repr = p._dict_repr
+        assert 'EasyBilayer' in dict_repr
+        assert 'front_head_layer' in dict_repr['EasyBilayer']
+        assert 'front_tail_layer' in dict_repr['EasyBilayer']
+        assert 'back_tail_layer' in dict_repr['EasyBilayer']
+        assert 'back_head_layer' in dict_repr['EasyBilayer']
+        assert 'constrain_heads' in dict_repr['EasyBilayer']
+        assert 'conformal_roughness' in dict_repr['EasyBilayer']
+
+
+def test_dict_round_trip():
+    """Test serialization/deserialization round trip."""
+    # When
+    d2o = Material(sld=6.36, isld=0, name='D2O')
+    air = Material(sld=0, isld=0, name='Air')
+
+    front_head = LayerAreaPerMolecule(
+        molecular_formula='C10H18NO8P',
+        thickness=12.0,
+        solvent=d2o,
+        solvent_fraction=0.3,
+        area_per_molecule=50.0,
+        roughness=2.0,
+        name='Custom Front Head',
+    )
+    tail = LayerAreaPerMolecule(
+        molecular_formula='C32D64',
+        thickness=18.0,
+        solvent=air,
+        solvent_fraction=0.0,
+        area_per_molecule=50.0,
+        roughness=2.0,
+        name='Custom Tail',
+    )
+    back_head = LayerAreaPerMolecule(
+        molecular_formula='C10H18NO8P',
+        thickness=12.0,
+        solvent=d2o,
+        solvent_fraction=0.4,
+        area_per_molecule=50.0,
+        roughness=2.0,
+        name='Custom Back Head',
+    )
+
+    p = Bilayer(
+        front_head_layer=front_head,
+        tail_layer=tail,
+        back_head_layer=back_head,
+        constrain_heads=False,
+        conformal_roughness=False,
+    )
+    p_dict = p.as_dict()
+    global_object.map._clear()
+
+    # Then
+    q = Bilayer.from_dict(p_dict)
+
+    # Expect
+    assert sorted(p.as_dict()) == sorted(q.as_dict())
+
+
+def test_dict_round_trip_constraints_enabled():
+    """Test round trip with constraints enabled."""
+    # When
+    p = Bilayer(constrain_heads=True, conformal_roughness=True)
+    p_dict = p.as_dict()
+    global_object.map._clear()
+
+    # Then
+    q = Bilayer.from_dict(p_dict)
+
+    # Expect
+    assert q.constrain_heads is True
+    assert q.conformal_roughness is True
+    assert sorted(p.as_dict()) == sorted(q.as_dict())
+
+
+def test_dict_round_trip_constraints_disabled():
+    """Test round trip with constraints disabled."""
+    # When
+    p = Bilayer(constrain_heads=False, conformal_roughness=False)
+    p_dict = p.as_dict()
+    global_object.map._clear()
+
+    # Then
+    q = Bilayer.from_dict(p_dict)
+
+    # Expect
+    assert q.constrain_heads is False
+    assert q.conformal_roughness is False
+    assert sorted(p.as_dict()) == sorted(q.as_dict())
+
+
+def test_constrain_multiple_contrast():
+    """Test multi-contrast constraints between bilayers."""
+    # When
+    p1 = Bilayer(name='Bilayer 1', constrain_heads=False)
+    p2 = Bilayer(name='Bilayer 2', constrain_heads=False)
+
+    # Set initial values
+    p1.front_head_layer.thickness.value = 10.0
+    p1.front_tail_layer.thickness.value = 16.0
+
+    # Constrain p2 to p1
+    p2.constrain_multiple_contrast(
+        p1,
+        front_head_thickness=True,
+        tail_thickness=True,
+    )
+
+    # Then - p2 values should match p1
+    assert p2.front_head_layer.thickness.value == 10.0
+    assert p2.front_tail_layer.thickness.value == 16.0
+
+    # Change p1 - p2 should follow
+    p1.front_head_layer.thickness.value = 12.0
+    assert p2.front_head_layer.thickness.value == 12.0
+
+
+def test_constrain_solvent_roughness():
+    """Test constraining solvent roughness to bilayer roughness."""
+    # When
+    p = Bilayer(conformal_roughness=True)
+    layer = Layer()
+
+    p.front_head_layer.roughness.value = 4.0
+
+    # Then
+    p.constrain_solvent_roughness(layer.roughness)
+
+    # Expect
+    assert layer.roughness.value == 4.0
+
+    # Change bilayer roughness - solvent should follow
+    p.front_head_layer.roughness.value = 5.0
+    assert layer.roughness.value == 5.0
+
+
+def test_constrain_solvent_roughness_error():
+    """Test error when constraining solvent roughness without conformal roughness."""
+    import pytest
+
+    p = Bilayer(conformal_roughness=False)
+    layer = Layer()
+
+    with pytest.raises(ValueError, match='Roughness must be conformal'):
+        p.constrain_solvent_roughness(layer.roughness)