def mutation_density(mut_ct_dict, neighbours):
    """
    Find the mutation density at each residue position

    Parameters
    ----------
    mut_ct_dict: dict
        dictionary containing full residue ids as keys
        and mutation counts as values
    neighbours: dict
        dictionary containing full residue ids as keys
        and lists of neigbours' residue ids as values

    Returns
    -------
    density_dict : dict
        Maps each residue to a number of mutations within
        its neighbors and itself. Neighbors are defined
        by an angstrom radius cutoff.
    """

    density_dict = {}
    for Residue1 in mut_ct_dict:
        # count muations on the residue itself
        num_mutations = mut_ct_dict[Residue1]
        # count mutations on neighbouring residues
        for Residue2 in neighbours[Residue1]:
            # will be present in mut_ct_dict
            # only if there was a mutation
            if Residue2 in mut_ct_dict:
                num_mutations += mut_ct_dict[Residue2]
        density_dict[Residue1] = num_mutations

    return density_dict