<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def define_density_matrix(Ne, explicitly_hermitian=False, normalized=False, variables=None):
 r"""Return a symbolic density matrix. The arguments are Ne (integer):
 The number of atomic states. explicitly_hermitian (boolean):
 Whether to make $\rho_{ij}=\bar{\rho}_{ij}$ for $i<j$ normalized (boolean):
 Whether to make $\rho_{11}=1-\sum_{i>1} \rho_{ii}$ A very simple example: Matrix([ [rho11, rho12], [rho21, rho22]]) The density matrix can be made explicitly hermitian Matrix([ [rho11, conjugate(rho21)], [rho21, rho22]]) or normalized Matrix([ [-rho22 + 1, rho12], [ rho21, rho22]]) or it can be made an explicit function of given variables Matrix([ [rho11(t, z), rho12(t, z)], [rho21(t, z), rho22(t, z)]]) """

    if Ne > 9:
        comma = ","
        name = r"\rho"
        open_brace = "_{"
        close_brace = "}"
    else:
        comma = ""
        name = "rho"
        open_brace = ""
        close_brace = ""

    rho = []
    for i in range(Ne):
        row_rho = []
        for j in range(Ne):
            if i == j:
                row_rho += [define_symbol(name, open_brace, comma, i, j,
                                          close_brace, variables,
                                          positive=True)]
            elif i > j:
                row_rho += [define_symbol(name, open_brace, comma, i, j,
                                          close_brace, variables)]
            else:
                if explicitly_hermitian:
                    row_rho += [conjugate(define_symbol(name, open_brace,
                                                        comma, j, i,
                                                        close_brace,
                                                        variables))]
                else:
                    row_rho += [define_symbol(name, open_brace, comma, i, j,
                                              close_brace, variables)]
        rho += [row_rho]

    if normalized:
        rho11 = 1-sum([rho[i][i] for i in range(1, Ne)])
        rho[0][0] = rho11

    rho = Matrix(rho)
    return rho