Input uncertainty | Installation | Quick Start | API | Chromatin Dynamics | Chromatin Mechanics

HIPPS-DIMES is a Python implementation of a maximum-entropy method for constructing ensembles of 3D chromatin structures from experimentally measured contact maps or pairwise distance constraints.¹²³ It accepts a mean distance map, a mean-squared-distance map, or a Hi-C contact map (converted internally into distance constraints) and generates conformations whose ensemble statistics fit those constraints under the selected optimization method.

HIPPS-DIMES can account explicitly for uncertainty in the input data. Its noise-aware maximum-entropy formulation weights agreement with each input target by its specified uncertainty, allowing noisy measurements to be fitted without requiring exact agreement. The COV method optimizes this formulation.

In addition to reconstructing static 3D chromatin structures, the model predicts chromatin-locus dynamics and chromatin-locus mechanics using polymer physics and the Ornstein–Uhlenbeck process. Available observables include autocorrelation functions (ACFs), mean-square displacements (MSDs), system-level viscoelastic moduli, and per-locus mechanical susceptibilities.

The theory and applications of this method are described in the following publications:

• Shi, Guang, and D. Thirumalai. "Epigenetic state encodes locus-specific chromatin mechanics." bioRxiv 2025.12.27.696709 (2025). link
• Shi, Guang, and D. Thirumalai. "From Hi-C contact map to three-dimensional organization of interphase human chromosomes." Physical Review X 11.1 (2021): 011051. link
• Shi, Guang, and D. Thirumalai. "A maximum-entropy model to predict 3D structural ensembles of chromatin from pairwise distances with applications to interphase chromosomes and structural variants." Nature Communications 14.1 (2023): 1150. link
• Shi, Guang, Sucheol Shin, and D. Thirumalai. "Static three-dimensional structures determine fast dynamics between distal loci pairs in interphase chromosomes." Science Advances 11.31 (2025): eadx1763. link

Other applications of this method can be found in:

• Dey, Atreya, et al. "Structural changes in chromosomes driven by multiple condensin motors during mitosis." Cell Reports 42.4 (2023).
• Jeong, Davin, et al. "Structural basis for the preservation of a subset of topologically associating domains in interphase chromosomes upon cohesin depletion." eLife 12 (2024): RP88564.

Release 4.0

HIPPS-DIMES 4.0 adds support for uncertainty in experimental data through a noise-aware maximum-entropy objective. The new COV method optimizes this objective on the CPU or GPU and checks convergence automatically.

This release also makes missing-data handling consistent across contact, distance, and squared-distance inputs. Missing pairs can be interpolated or excluded, and fully missing loci require an explicit repair or removal policy.

Accounting for uncertainty in input data

Experimental contact maps and pairwise distance measurements contain noise. Treating their values as exact constraints can force a model to fit measurement errors; the targets may even be incompatible with any physical ensemble. Noise-aware HIPPS-DIMES retains the maximum-entropy principle while allowing deviations from the targets according to their specified uncertainties. More precise observations constrain the fit more strongly, while less precise observations permit larger deviations.

Bayesian formulation and objective

Let $D_{ij}^{\mathrm{obs}}$ denote a mean-squared-distance target and let $D_{ij}(B)=B_{ii}+B_{jj}-2B_{ij}$ be its model prediction. Here, $B_{ij}=\langle\mathbf r_i\cdot\mathbf r_j\rangle$ is the ensemble-averaged Gram matrix of coordinates measured relative to the center of mass. We assume independent Gaussian measurement errors with specified variances $v_{ij}>0$:

$$ D_{ij}^{\mathrm{obs}}=D_{ij}(B)+\varepsilon_{ij}, \qquad \varepsilon_{ij}\sim\mathcal N(0,v_{ij}). $$

To retain the preference for maximum entropy, we choose an entropy-favoring prior $\pi(B)\propto\exp[S(B)]$ over physically valid Gram matrices, where $S(B)$ is the Gaussian ensemble entropy in units of Boltzmann's constant. Bayes' theorem combines this prior with the Gaussian measurement likelihood:

$$ p(B\mid D^{\mathrm{obs}}) \propto p(D^{\mathrm{obs}}\mid B),\pi(B). $$

With the variances held fixed, maximizing this posterior is equivalent to minimizing its negative logarithm. Since $S(B)=\frac32\log\det{}'B+\mathrm{constant}$ in three dimensions, the resulting objective, up to an additive constant, is

$$ F(B)=-\frac32\log\det{}'B +\frac12\sum_{(i,j)\in\mathcal O} \frac{\left[D_{ij}(B)-D_{ij}^{\mathrm{obs}}\right]^2}{v_{ij}}. $$

The first term favors higher ensemble entropy. The second penalizes disagreement with the data, weighted by each pair's uncertainty. The set $\mathcal O$ contains each included pair once ($i<j$). The determinant $\det{}'B$ is over the $N-1$ internal modes: overall translation is removed, and the internal eigenvalues must remain positive.

The noise model applies to the final mean-squared-distance targets, after input conversion, normalization, and missing-data handling. Interpolated or repaired pairs enter the fit with the chosen noise model; explicitly excluded pairs do not. The user supplies the noise level, which is held fixed during optimization.

Optimizing the objective with COV

COV minimizes this objective by optimizing the centered Gram matrix $B$ while keeping its internal eigenvalues positive. Its default hybrid optimizer combines primal-dual hybrid gradient (PDHG) with the fast iterative shrinkage-thresholding algorithm (FISTA). Both stages optimize the same objective, and an independent optimality check determines convergence. CPU and float64 GPU execution are supported.

This gives a maximum a posteriori (MAP) estimate under the stated entropy-favoring prior and Gaussian error model. COV returns one fitted model; generated conformations sample that model's structural ensemble, not a posterior over model parameters.

See the COV quick start for usage, noise-model notes for the uncertainty conventions and initialization, and solver documentation for algorithmic details.

Documentation

Requirements

• Python 3.9+

Install

First, clone this repository:

git clone https://github.com/anyuzx/HIPPS-DIMES
cd HIPPS-DIMES

Using pip

pip install -e .

Using uv

uv pip install -e .

Either command installs the dependencies declared in pyproject.toml, installs HIPPS-DIMES as a Python package, and provides the HippsDimes command-line entry point.

Dependencies

The package requires:

• Python 3.9+
• Click - Command-line interface
• rich-click - Rich-styled command-line interface
• NumPy - Numerical computing
• SciPy - Scientific computing
• pandas - Data manipulation
• tqdm - Progress bars
• cooler - .mcool and .cool data-format support
• Rich - Rich terminal output
• hic-straw - .hic data-format support

Optional (for GPU acceleration):

• CuPy - GPU-accelerated computing via CUDA

To install CuPy for GPU support:

conda install -c conda-forge cupy
# or
pip install cupy-cuda11x  # Replace with your CUDA version

Quick start

The quickest way to get started is with one of the Google Colab notebooks:

• Quick start
• Chromatin Dynamics and Mechanics

To view the available arguments and options, run:

HippsDimes --help

COV quick start

Use --method COV to account for input uncertainty. Specify exactly one noise model:

• --gaussian-noise-variance <v> assigns the same absolute variance to every included pair, $v_{ij}=v$ (homoskedastic uncertainty).
• --gaussian-noise-relative-std <r> assigns the same relative standard deviation to every included pair, $\sigma_{ij}=rD_{ij}^{\mathrm{obs}}$. The absolute variances then depend on the pair, $v_{ij}=(rD_{ij}^{\mathrm{obs}})^2$ (heteroskedastic uncertainty).

For example, fit a mean-squared-distance map with an assumed relative standard deviation of 10%:

python -m hipps_dimes observed_ddmap.npy cov_fit \
  --input-type ddmap --input-format npy \
  --method COV --gaussian-noise-relative-std 0.1 \
  --iteration 10000

Here, 0.1 means that the standard deviation is 10% of each final mean-squared-distance target. It is an illustrative assumption supplied by the user, not an estimated noise level or 10% uncertainty in the original contact counts. The command runs on the CPU; add --use-gpu if CuPy and an accessible CUDA GPU are available.

COV stops when its convergence test passes or the iteration budget is exhausted. Check the reported convergence status before using the fit; see COV convergence and returned results. For incomplete data, see the missing-data policies and observation requirements.

API

Input files

HIPPS-DIMES accepts contact maps (cmap), mean distance maps (dmap), and mean-squared-distance maps (ddmap). Contact maps may use .cool or .mcool format through cooler, .hic format through hicstraw, a plain-text matrix, or a NumPy .npy matrix. Mean distance and mean-squared-distance maps may use plain text or NumPy .npy format.

In a plain-text matrix, each line contains one space-separated row. For example, a mean distance map could look like:

0  2  3
2  0  2
3  2  0

Output files

This script will generate several files:

• A text file containing the final model mean distance map.
• For contact-map inputs, a text file containing the model contact map at the optimized contact threshold. This threshold maximizes agreement with the normalized experimental contact map.
• For cooler and .hic contact-map inputs, a text file containing the internal target contact map used to construct the distance constraints: {output_prefix}_cmap_target.txt.
• A text file containing the connectivity matrix.
• An .xyz file containing the generated ensemble (optional).
• A CSV file containing run parameters: {output_prefix}_run_parameters.csv (optional).
• A CSV file containing iteration-series scalar data: {output_prefix}_iteration_series.csv (optional).

Explanation of the arguments and options

Arguments

• INPUT: Path to a contact map, mean distance map, or mean-squared-distance map.
• OUTPUT_PREFIX: Prefix for output files. For example, specifying TEST makes the output filenames start with TEST_.

Options

• -k, --connectivity-matrix: Path to a connectivity matrix to use as the initialization, for example when restarting from a previous result.
• -e, --ensemble: Number of individual conformations to generate from the fitted model. Default: 1000.
• -a, --alpha: Contact-to-distance conversion exponent. For contact-map input, HIPPS-DIMES uses $d_{ij} \propto c_{ij}^{-1/\alpha}$. The default value, $\alpha=4.0$, was estimated in this work. Default: 4.0.
• -s, --selection: Specify chromosome or region. This option is required when the input file has cooler or .hic format. For cooler files, the value is passed to the cooler.Cooler.matrix().fetch() method. For .hic files, use format "chr1:start1-end1,chr2:start2-end2". For details on cooler selectors, please refer to their documentation.
• -m, --method: Select IS (Iterative Scaling, default), GD (Gradient Descent), DI (Direct Inversion), or COV (optimization of the noise-aware maximum-entropy objective).
• -l, --lamd: L1 or L2 regularization weight. A value of 0.0 disables regularization. Must remain 0.0 for COV, whose noise model determines the uncertainty weighting. Default: 0.0.
• -r, --reg: Regularization type: L1 or L2 (default). Use this option together with --lamd.
• --gaussian-noise-variance: Positive scalar absolute variance shared by all included mean-squared-distance targets. COV only.
• --gaussian-noise-relative-std: Positive scalar relative standard deviation sigma_ij / Dobs_ij. COV converts it to pair variance (value * Dobs_ij)^2, using the final mean-squared-distance targets after preprocessing and input conversion. COV requires exactly one of this option and --gaussian-noise-variance.
• --covariance-optimizer: COV optimizer: hybrid (default) or pdhg. The hybrid combines PDHG with FISTA refinement; pdhg uses PDHG alone. See the solver documentation for advanced tuning.
• --covariance-relative-tolerance: Relative tolerance for the COV optimality (KKT) test. Default: 1e-5.
• --covariance-absolute-tolerance: Absolute tolerance for internal COV optimality checks. Default: 1e-10.
• --covariance-handoff-relative-tolerance: Relative optimality threshold for switching from PDHG to FISTA in the hybrid optimizer. Default: 1e-2.
• -i, --iteration: Maximum optimizer iterations. Default: 10000.
• --learning-rate: Learning rate for IS or GD. Typical IS values are 1–30; GD generally requires a much smaller value, such as 1e-8. Default: 10.0.
• --momentum: Momentum coefficient for IS, between 0.0 and 1.0. Recommended: use 0.95 with --nesterov for the fastest convergence observed in benchmarks. Use 0.9 for a more conservative setting. Default: 0.0.
• --nesterov: Use Nesterov Accelerated Gradient with IS. Recommended with --momentum 0.95.
• --use-gpu: Enable GPU acceleration through CuPy. All COV optimizers use float64. COV fails rather than silently falling back when CUDA is unavailable.
• --gpu-float32: Use float32 for legacy GPU IS/GD. COV is float64-only.
• --save-steps: Comma-separated list of iteration steps at which to save the connectivity matrix. Example: --save-steps 1000,5000,10000. Files are saved as {output_prefix}_connectivity_matrix_iter{step}.txt. When used as a library (without output_prefix), connectivity matrices at these steps are still returned in results['connectivity_matrix_at_steps'].
• --eigh-threads: Number of eigendecomposition and BLAS/LAPACK threads. If unset, the backend default is used. Set to 1 for single-threaded runs.
• --input-type: Required input type: cmap (contact map), dmap (mean distance map), or ddmap (mean-squared-distance map).
• --input-format: Required input format: text, npy, cooler, or hic. Contact maps support all four formats; dmap and ddmap support text and npy.
• --binsize: Bin size for .hic input, in bp. Default: 25000.
• --norm: .hic normalization: KR, VC, or NONE. Default: KR.
• --unit: .hic unit: BP or FRAG. Default: BP.
• --no-log: By default, the program writes two log files when output-prefix is provided: {output_prefix}_run_parameters.csv and {output_prefix}_iteration_series.csv. Use --no-log to disable writing both files.
• --no-xyzs: Do not write generated conformations to an .xyz file.
• --ignore-missing-data: Exclude remaining missing pair constraints. Without this flag, every remaining missing pair is interpolated before optimization: contact maps in log-contact space, distance maps in distance space, and squared-distance maps in squared-distance space. The completed map is the optimization target.
• --repair-fully-missing-loci: Impute only the genomic nearest-neighbor constraints (i, i-1) and (i, i+1) needed to reconnect a locus with no observed off-diagonal pairs. The repaired values become ordinary target constraints in the selected objective and, for COV, in its noise model.
• --remove-fully-missing-loci: Remove loci that have no observed off-diagonal pairs before applying the remaining missing-pair policy.
• --balance: Balance a cooler-format contact map before optimization.
• --neighbor-balance: Apply neighbor balancing to a contact map by dividing each pair value by the geometric mean of the corresponding neighbor-contact values. See Paggi and Zhang (2025) for details.
• --not-normalize: Disable automatic maximum-value normalization of a contact map.
• --enforce-nonnegative-connectivity-matrix: Constrain all off-diagonal spring constants to be nonnegative. Cannot be combined with COV.
• --save-pickle: Save the returned results dictionary as {output_prefix}_HIPPS_DIMES_results.pkl instead of writing the default text, CSV, and XYZ outputs.
• -q, --quiet: Disable table output while retaining the progress bar.

Examples

The following input-format examples use the default IS method, which does not model input uncertainty. For noise-aware fitting, use COV with one of the noise options shown in the COV quick start.

Example 1

First, download a cooler-format Hi-C contact map from here (the file size is about 116 MB). This map represents a mitotic chicken chromosome and was originally retrieved from the GEO repository. Rename it hic_example.cool, and then run:

HippsDimes hic_example.cool test --input-type cmap --input-format cooler -s chr7:10M-15M -i 10 -e 10

This command loads hic_example.cool and runs Iterative Scaling. The test argument makes output filenames start with test_; --input-type cmap identifies the input as a contact map; --input-format cooler identifies the file format; and -s chr7:10M-15M selects the 10–15 Mb region of chromosome 7. The input type, input format, and cooler region selection are required here. Run HippsDimes --help for the complete option requirements.

With these defaults, the program writes test.xyz, test_connectivity_matrix.txt, test_dmap_final.txt, test_cmap_final.txt, test_cmap_target.txt, test_run_parameters.csv, and test_iteration_series.csv. The test.xyz file contains 10 conformations and can be viewed with VMD or other compatible visualization software.

Example 2

This example uses a Hi-C contact map for chromosome 14 in HeLa cells, measured 12 hours after release from prometaphase. Download the .cool file from here which was originally retrieved from the GEO repository under accession GSE102740. The file size is about 655 MB. Once downloaded, run:

HippsDimes GSM3909682_TB-HiC-Dpn-R2-T12_hg19.1000.multires.cool::6 test --input-type cmap --input-format cooler -s chr14:20M-107M -i 10000 -e 10

This command loads group 6 of the multiresolution cooler file and runs HIPPS-DIMES for 10,000 iterations. On an AMD Ryzen 5 3600 CPU, this example takes approximately 3–4 minutes.

Example 3

A .hic example:

HippsDimes mydata.hic test \
  --input-type cmap --input-format hic \
  --selection chr1:31000000-41000000,chr1:31000000-41000000 \
  --binsize 25000 --norm KR --unit BP \
  -i 10000 -e 10

For an example of applying HIPPS-DIMES directly to imaging data, see the corresponding notebook.

Tips

• Computational cost and memory use grow rapidly with matrix size, but there is no universal matrix-size cutoff for convergence. The practical limit depends on the selected optimizer, hardware, number of observed pairs, and iteration budget. If a problem is too expensive, coarse-grain the input or analyze a smaller genomic region.
• Optimization tuning is method-specific. --learning-rate applies to IS and GD, while --momentum and --nesterov apply only to IS. Typical IS learning rates are between 1 and 30; GD generally requires a much smaller value, such as 1e-8. These parameters are not used by DI or COV. Nesterov acceleration, for example --momentum 0.95 --nesterov, can accelerate IS, but the benefit is problem-dependent.
• For larger problems, consider --use-gpu when CuPy and an accessible CUDA GPU are available. The speedup depends on the matrix, optimizer, GPU, and CPU, so benchmark a representative case rather than assuming a fixed acceleration factor. COV uses float64 on the GPU and does not silently fall back to the CPU.
• Use --save-steps 1000,5000,10000 to retain intermediate connectivity matrices. For IS, these checkpoints support the manual entropy-and-loss stopping decision described below. For COV, the built-in KKT test determines convergence; checkpoints are useful for diagnostics or restarting but do not replace the reported convergence certificate.
• Missing pairs are interpolated by default and the completed map becomes the target. Use --ignore-missing-data to exclude them instead. Non-finite input pairs are missing for every map type; nonpositive off-diagonal pairs are also missing for contact maps. If a locus has no observed off-diagonal pairs at all, select either --repair-fully-missing-loci to add its genomic-neighbor constraints or --remove-fully-missing-loci to remove it before the pair policy is applied, regardless of --ignore-missing-data. Repair and removal are mutually exclusive. For COV, the included pairs must connect all retained loci, so disconnected groups are rejected even when no individual locus is isolated. Their relative motion would otherwise be unconstrained, leaving the objective without a finite optimum. DI requires a complete target and cannot be combined with excluded pairs.
• Contact-map inputs are normalized by their maximum entry by default. Use --not-normalize when the supplied contact map should retain its existing scale.
• A contact map alone does not define a physical length scale, so its inferred distances and structures are dimensionless. An external distance measurement, such as the mean distance between neighboring loci, can be used to rescale them.

How to choose the optimization method, iteration count, and convergence criterion

The recommended method and stopping rule depend on the input and on whether its uncertainty should be modeled:

Input and interpretation	Recommended method	How to stop
Experimental contact map, with uncertainty	COV with a user-specified noise model and noise level	Automatic KKT convergence test; --iteration is the maximum update budget
Experimental contact map, uncertainty not modeled	IS is a practical alternative when only static 3D structures are needed	Manual judgment from the entropy and loss histories
Valid mean-squared-distance map, treated as exact	DI (preferred) or IS	DI is a one-step conversion; IS requires enough iterations to reproduce the target
Valid mean-squared-distance map, with uncertainty	COV with a user-specified noise model and noise level	Automatic KKT convergence test; --iteration is the maximum update budget

Experimental contact maps

For experimental contact maps, the preferred method is COV with a noise model and uncertainty level supplied by the user. It fits the converted mean-squared-distance targets according to their uncertainties using the noise-aware objective. See the COV quick start for the two supported noise models and convergence guidance for checking a fit.

Converting an experimental contact map to pairwise distances does not guarantee that the resulting mean-squared-distance map is a valid Euclidean distance matrix. Consequently, the exact hard-constraint problem approached by Iterative Scaling (IS) may have no feasible global solution: no Gaussian ensemble may satisfy all inferred pair distances simultaneously.

When uncertainty is not modeled and only static 3D structures are needed, IS remains a practical alternative. IS has no built-in convergence criterion and always runs for the number of iterations requested by the user. The iteration-series file records both entropy and loss; here, loss is the root-mean-square relative difference between the model and target mean-squared distances over the constrained pairs. As a default, select the iteration with the highest recorded entropy. In the rarer case where the entropy curve has an early or local maximum while the loss is still changing materially, do not treat that peak as convergence. Continue until the loss begins to level off, and select a checkpoint using both curves. This is necessarily an educated user judgment rather than a convergence certificate supplied by HIPPS-DIMES.

Warning — IS with experimental contact maps: Because the inferred target may have no feasible global solution, the high-frequency relaxation modes continue to drift as more IS iterations are performed. Interpret these modes with caution. This sensitivity is most relevant to the high-frequency regime of the predicted loss modulus, $G''(\omega)$. The low-frequency modes converge and, once converged, do not depend on the selected number of IS iterations.

Valid mean-squared-distance maps

A mean-squared-distance map computed directly from an ensemble of 3D coordinates is a valid Euclidean distance matrix, up to numerical precision. If it is to be treated as exact, use Direct Inversion (--method DI). DI performs the covariance-to-connectivity conversion in one step, so no iteration count or iterative convergence decision is needed.

For a complete valid target, IS with suitable numerical settings converges to the same solution as DI, but the required number of iterations depends on the system size and optimization settings. DI is therefore preferred when the map is valid and uncertainty is intentionally ignored. If uncertainty in the mean-squared distances should be represented, use COV instead and specify either an absolute variance or a relative standard deviation.

COV convergence and returned results

COV checks whether the fitted model satisfies the objective's first-order optimality conditions, also called the Karush–Kuhn–Tucker (KKT) conditions. Before reporting convergence, it independently recomputes this check on the returned model. The default relative tolerance is 1e-5.

--iteration sets the maximum number of updates. Reaching that limit does not establish convergence. If the returned model fails the convergence check, the command-line program retains requested output files but exits with status 1. The Python API emits a RuntimeWarning and returns the partial result with results["covariance_optimization"]["converged"] == False.

In the Python API, inspect results["covariance_optimization"]:

• converged reports whether the independent check passed.
• status explains why the optimizer stopped, and relative_eliminated_kkt_residual gives the relative optimality residual.
• iterations counts executed updates; returned_iteration identifies the selected model, which may come from an earlier update.

The iteration-series column is_returned_iterate marks the selected model. The reported COV objective, loss, and entropy all refer to that model. The optimized objective combines entropy and uncertainty-weighted data fit; loss is a separate diagnostic of relative distance error. See the solver documentation for detailed diagnostics and tuning.

Use HIPPS-DIMES as a Python library

HIPPS-DIMES can be used as either a command-line tool or a Python library, making it straightforward to integrate into Python workflows, Jupyter notebooks, and automated pipelines.

Main function: run_optimization()

The core functionality is available through the run_optimization() function. This example uses COV with an illustrative 10% relative standard deviation on the converted mean-squared-distance targets. Choose a noise level appropriate to your data.

import numpy as np
import hipps_dimes as HD

# Load your contact map
cmap = np.loadtxt("contact_map.txt")

# Fit while accounting for input uncertainty
results = HD.run_optimization(
    input_matrix=cmap,          # Provide matrix directly
    input_type="cmap",
    method="COV",
    gaussian_noise_relative_std=0.1,
    iteration=10000,
    ensemble=1000,
    verbose=False,              # Suppress console output
)

if not results["covariance_optimization"]["converged"]:
    raise RuntimeError("COV fit did not converge; inspect its diagnostics.")

# Access results
connectivity_matrix = results["connectivity_matrix"]
structures = results["xyzs"]          # (ensemble, n_beads, 3)
final_dmap = results["dmap_final"]
final_cmap = results["cmap_final"]
iteration_series = results["iteration_series"]
run_parameters = results["run_parameters"]
optimization = results["optimization"]

Add use_gpu=True to run on an accessible CUDA GPU with CuPy installed. Omitting method still selects IS; COV must be requested explicitly.

Return values

The run_optimization() function returns a dictionary with:

• 'connectivity_matrix': Final connectivity matrix (NumPy array)
• 'dmap_final': Final distance map (NumPy array)
• 'cmap_final': Final contact map (NumPy array, for input_type="cmap")
• 'xyzs': Generated conformations (NumPy array, unless no_xyzs=True)
• 'iteration_series': Iteration-series scalar outputs (pandas DataFrame). Every method reports iteration, loss, and entropy; each numbered row describes the model produced by that solver step. In particular, IS/GD row t is the state after update t. COV adds optimizer-specific diagnostics.
• 'run_parameters': Run parameters (pandas DataFrame with parameter and value columns)
• 'optimization': Common returned-model summary for every method, containing method, status, converged, iterations_executed, returned_iteration, final_loss, and final_entropy. IS and GD report converged=None because they run a fixed iteration budget without a convergence certificate; DI reports converged=None because it is a direct solve. COV reports its independently certified Boolean convergence status.
• 'log': Alias for 'iteration_series' (backward compatibility)
• 'rc_optimal': Optimal contact threshold (float, for input_type="cmap")
• 'connectivity_matrix_at_steps': Saved intermediate connectivity matrices (dictionary, when save_steps is set)
• 'gram_matrix' and 'covariance_optimization': Fitted Gram matrix and convergence diagnostics (for method="COV")

Additional utility functions

The package also provides helper functions for direct use:

import hipps_dimes as HD

# Generate structures from connectivity matrix
structures = HD.a2xyz_sample(connectivity_matrix, ensemble=1000)

# Compute an iteratively aligned mean structure
mean_structure, iterations, difference = HD.iterative_generalized_procrustes_alignment(
    structures, tolerance=1e-8, max_iterations=100, allow_reflection=True
)

# Convert connectivity matrix to distance map
dmap = HD.a2dmap_theory(connectivity_matrix)

# Convert connectivity matrix to contact map
cmap = HD.a2cmap_theory(connectivity_matrix, rc=5.0)

# Create a Rouse chain connectivity matrix
A = HD.construct_connectivity_matrix_rouse(n=100, k=1.0)

iterative_generalized_procrustes_alignment accepts coordinates shaped (n_structures, n_points, n_dimensions) with corresponding points in the same order. It centers each structure, initializes the reference with the first structure, and repeatedly uses optimal_rotate to align the original centered structures to the current mean. Reflections are allowed by default; set allow_reflection=False for proper rotations only. No scaling is applied, and the input is not modified. It returns the centered mean, iteration count, and final RMS change over all point coordinates. Convergence means difference < tolerance; if the iteration limit is reached first, the latest mean and change are returned.

---

Dynamics Prediction Functionality

In addition to reconstructing static 3D chromatin structures from contact or distance maps, HIPPS-DIMES can simulate chromatin dynamics.³ The dynamics are based on polymer physics and the Ornstein–Uhlenbeck process and provide time-dependent observables such as autocorrelation functions (ACFs) and the mean-square displacements (MSDs) of individual loci.

Functions

• compute_acf_general_theory(i, j, t, a, zeta=1.0): Numerically computes the time-dependent autocorrelation function between monomers i and j from a connectivity matrix a. It also returns the corresponding two-point MSD for every time in t.

• compute_m1_i(i, t, a, zeta=1.0): Computes the single-locus MSD for monomer i. The returned two-dimensional array contains time in the first column and MSD in the second.

Dynamics class

The Dynamics class provides trajectory simulation from a connectivity matrix a.

Example code

import hipps_dimes as HD

model = HD.Dynamics(a)  # a is the connectivity matrix
model.initialize(dt=1e-2, zeta=1.0, beta=1.0)

model.run(int(1e5), every=10)
model.resume(int(5e4), every=10)

Trajectory coordinates are available in model.traj, a (T, N, 3) NumPy array, where T is the number of snapshots and N is the number of loci. The reduced simulation time for each saved snapshot is stored in model.traj_time, a length-T NumPy array.

Save both arrays with model.save_traj("traj.npz"). The .npz file contains the traj and traj_time arrays.

Dynamics.run(...) starts a fresh trajectory and can be called only once per simulation state. To continue an existing simulation, use Dynamics.resume(...). When arguments are omitted, resume(...) reuses the previous passive simulation settings for update, every, method, and update_zero_modes; any of them may still be overridden explicitly. To discard the previous trajectory and start over on the same object, call Dynamics.reset() before run(...). By default, run(...) does not append the post-integration final state to model.traj; set include_final_state=True to include it.

Dynamics under external force: Dynamics.run_with_force(...)

In addition to passive dynamics (Dynamics.run), you can simulate trajectories with a constant external force applied to selected loci.

Key parameters

• force_loci: list of locus indices where the force is applied
• force_amplitude: force magnitude
• force_direction: (3,) direction vector (normalized internally)
• force_duration: optional number of time steps to apply the force; if None, the force is applied for the entire run

Example: forced dynamics

import numpy as np
import hipps_dimes as HD

# Load a connectivity matrix or obtain one from HD.run_optimization().
a = np.loadtxt("my_connectivity_matrix.txt")
model = HD.Dynamics(a)
model.initialize(dt=1e-2, zeta=1.0, beta=1.0)

# Pull locus 10 along +x for the first 2e4 steps (then release)
model.run_with_force(
    T=int(1e5),
    force_loci=[10],
    force_amplitude=1.0,
    force_direction=[1.0, 0.0, 0.0],
    force_duration=int(2e4),
    every=10,
)

traj = model.traj  # shape: (n_snapshots, N, 3)

Linear mechanical response

HIPPS-DIMES provides utilities to compute system-level linear viscoelastic moduli and per-locus mechanical susceptibilities from a connectivity matrix a. These routines decompose the polymer into normal modes and exclude the zero, or center-of-mass, mode.

Note on units: freq is interpreted as angular frequency $\omega$. The returned response functions are in the model's internal units and depend on the friction coefficient zeta used to define relaxation times.

compute_modulus(a, freq, zeta=1.0)

Computes system-level moduli by summing contributions from all nonzero normal modes.

• Inputs
  • a: (N, N) symmetric connectivity matrix
  • freq: (n_freq,) array of angular frequencies $\omega$
  • zeta: friction coefficient (default 1.0)
• Returns
  • G_storage: (n_freq, 2) array with columns $[\omega, G'(\omega)]$
  • G_loss: (n_freq, 2) array with columns $[\omega, G''(\omega)]$

compute_monomer_mechanical_susceptibility(a, freq, zeta=1.0)

Computes the real and imaginary parts of the per-locus mechanical susceptibility:

$$ \chi_i'(\omega) = \frac{1}{\zeta}\sum_{p>0}v_{pi}^2 \frac{\tau_p}{1+(\omega\tau_p)^2}, \qquad \chi_i''(\omega) = \frac{1}{\zeta}\sum_{p>0}v_{pi}^2 \frac{\omega\tau_p^2}{1+(\omega\tau_p)^2}. $$

Here, $\tau_p=-\zeta/\lambda_p$. There is no factor of two in these monomer-level response functions.

• Returns
  • freq: (n_freq,)
  • chi_prime_i: (n_freq, N) array containing $\chi_i'(\omega)$
  • chi_double_prime_i: (n_freq, N) array containing $\chi_i''(\omega)$

Example: compute mechanical response

This example uses a converged COV fit. As above, the 10% relative noise level is illustrative and applies to the converted mean-squared-distance targets.

import numpy as np
import hipps_dimes as HD

# Example: obtain a connectivity matrix 'a' from HIPPS-DIMES
# (you can also load a saved matrix from disk with np.loadtxt)
results = HD.run_optimization(
    input_matrix=np.loadtxt("contact_map.txt"),
    input_type="cmap",
    method="COV",
    gaussian_noise_relative_std=0.1,
    iteration=10000,
    verbose=False,
)
if not results["covariance_optimization"]["converged"]:
    raise RuntimeError("COV fit did not converge; inspect its diagnostics.")
a = results["connectivity_matrix"]

# (1) Compute bulk moduli G'(ω), G''(ω)
freq = np.logspace(-3, 3, 200)  # angular frequencies ω
G_storage, G_loss = HD.compute_modulus(a, freq, zeta=1.0)

# (2) Compute per-locus mechanical susceptibilities
freq_out, chi_prime_i, chi_double_prime_i = (
    HD.compute_monomer_mechanical_susceptibility(a, freq, zeta=1.0)
)

How to cite

If you use this program in a publication, please cite the following references:

• Shi, Guang, and D. Thirumalai. "From Hi-C Contact Map to Three-dimensional Organization of Interphase Human Chromosomes." Physical Review X 11.1 (2021): 011051.

• Shi, G., Thirumalai, D. A maximum-entropy model to predict 3D structural ensembles of chromatin from pairwise distances with applications to interphase chromosomes and structural variants. Nat Commun 14, 1150 (2023).

• Shi, G., Shin, S., and Thirumalai, D. "Static three-dimensional structures determine fast dynamics between distal loci pairs in interphase chromosomes." Science Advances 11.31 (2025): eadx1763.

Footnotes

1. Shi, Guang, and D. Thirumalai. "From Hi-C Contact Map to Three-dimensional Organization of Interphase Human Chromosomes." Physical Review X 11.1 (2021): 011051. ↩

2. Shi, G., Thirumalai, D. A maximum-entropy model to predict 3D structural ensembles of chromatin from pairwise distances with applications to interphase chromosomes and structural variants. Nat Commun 14, 1150 (2023). ↩

3. Shi, G., Shin, S., and Thirumalai, D. "Static three-dimensional structures determine fast dynamics between distal loci pairs in interphase chromosomes." Science Advances 11.31 (2025): eadx1763. ↩ ↩²